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ABSTRACT 


This paper describes a new interacting boundary- layer approach for computing the 
viscous transonic flow over airfoils. The theory includes a complete treatment of 
viscous-interaction effects induced by the wake and accounts for normal pressure 
gradient effects across the boundary layer near trailing edges. The method is based 
on systematic expansions of the full Reynolds equation of turbulent flow in the limit 
of Reynolds numbers, Re <*>. The theory employs a local inner solution to describe 
the strong interaction region at trailing edges. Procedures are developed for 
incorporating the local trailing-edge solution into the numerical solution of the 
coupled full-potential and integral boundary- layer equations. Although the theory is 
strictly applicable to airfoils with cusped or nearly-cusped trailing edges and to 
turbulent boundary layers that remain fully attached to the airfoil surface, the 
method has been successfully applied to more general airfoils and to flows with small 
separation zones. Comparisons of theoretical solutions with wind-tunnel data indi- 
cate the present method can accurately predict the section characteristics of 
airfoils including the absolute levels of drag. The results of the study clearly 
demonstrate the importance of including both wake and normal pressure-gradient 
effects in the theoretical formulation and of using a fully conservative difference 
scheme in the solution of the inviscid equations. This study indicates that a simple 
interacting boundary- layer approach can be an effective tool for the prediction of 
viscous effects on airfoils. 
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1. INTRODUCTION 


The development of theoretical methods for the prediction of the section charac- 
teristics of airfoils has always been an important problem area in aerodynamics. It 
continues to receive much attention today, with interest focusing on the numerical 
solution of various theoretical models for the viscous flow over airfoils. Recent 
advances in numerical techniques for solving the full non-linear potential equation 
have led to practical methods for computing the lift and drag of inviscid transonic 
flows containing shock waves (Refs. 1, 2 & 3). It is, however, well known that 
viscous effects are important and must be taken into account in these flows if accu- 
rate section characteristics are to be predicted. This is particularly true for 
supercritical flow over rear-loaded airfoils for which the combination of shock waves 
and aft camber combine to produce significantly thickened boundary layers over the 
rear upper surface. This in turn leads to much larger viscous effects than are expe- 
rienced on conventional airfoils in supercritical flow. For heavy rear loading, the 
boundary layer can reduce the lift by a factor of two below inviscid levels, even at 
Reynolds numbers as large as 10^, 

In most practical problems the Reynolds number is large, the boundary layers are 
mostly turbulent, and the direct effects of viscosity and turbulent transport are 
confined to thin shear layers on the airfoil surface and along the wake. In these 
situations the viscous flow can be effectively analyzed by interacting boundary- layer 
theory (IBLT) in which the flow field is divided into a primary inviscid region, thin 
shear layers, and localized strong-interaction regions. Previous applications of the 
theory to airfoils used a highly simplified formulation based on the neglect of pres- 
sure variations across the shear layers and neglect of all interaction effects 
induced by the wake. This formulation led to a description in terms of the coupled 
inviscid and Prandtl boundary- layer equations which are solved by standard iterative 
schemes. These theories account for the primary viscous effect due to the displace- 
ment thickness on the airfoil but do not properly treat viscous effects due to the 
wake and in shock-wave/boundary- layer and trailing-edge interaction regions. 
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The new formulation developed in the present study extends the conventional 
boundary- layer approach to include a full treatment of the wake and normal pres- 
sure-gradient effects at the trailing edge. Our method is based on systematic 
expansions of the full Reynolds equation of turbulent flow in the limit of Reynolds 
number, Re «>. The method of matched asymptotic expansions leads to a description 
of the flow in terms of inviscid regions and thin shear layers near the airfoil and 
wake. A local inner solution is developed to correct the standard Prandtl 
boundary- layer formulation for strong- interaction effects at the trailing edge. This 
work can be viewed as an extension and further development of Melnik and Chow's 
turbulent-interaction theory previously described in Refs. 4 & 5. The method is 
applicable for high Reynolds number flow over general airfoil shapes with cusped or 
nearly cusped trailing edges and free-stream Mach numbers less than one. The method 
is intended for flows that are turbulent over most of the airfoil and that are not 
separated. However, provisions were made to allow for the presence of small sepa- 
ration zones so that the resulting code will function and provide at least a rough 
description of the solution for these cases. 

The formulation of the present paper employs a version of the mixed-flow relax- 
ation technique developed by Jameson (Refs. 6 & 7) for solving the full-potential 
equation in conservation form*. The method employs a conformal mapping of the 
airfoil to a circle to obtain a useful computational grid and uses a direct solver to 
accelerate convergence. The inviscid boundary conditions are modified to account for 
viscous effects using a surface-source formulation of the matching conditions thus 
avoiding the need to carry out repeated conformal mappings as required in displace- 
ment-thickness approaches. An iterative scheme is employed to obtain a 
self-consistent solution of the coupled boundary- layer and inviscid-f low equations. 
The viscous matching conditions employed in the theory account for displacement 
effects on the airfoil as well as both wake-thickness and wake- curvature effects. 
The local inner solution developed in Ref. 5 is used to correct the standard 
boundary- layer solution for strong-interaction effects at the trailing edge. The 
modified viscous method accounts for the significant pressure variation across the 


*A nonconservative option is also provided in the computer code. 
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boundary layer and removes the nonuniformity of classical boundary- layer theory near 
trailing edges. The resulting theory provides a completely rational treatment of the 
trailing-edge region for cusped airfoils at lift. It does not account for the addi- 
tional nonuniformity arising on airfoils with sharp but non-cusped trailing edges. 
However, this nonuniformity can be ignored and the resulting computer code can be 
applied to general airfoils with nonzero included angles. Results given later in the 
paper indicate that the method gives good overall results for airfoils with small 
trailing-edge angles (less than 10°). 

The present theory also does not provide for a proper treatment of the 
strong-interaction phenomenon near shock-wave/boundary- layer interaction zones. It 
ignores the fact that the boundary- layer approximations fail in these regions and 
determines the solution with a conventional interacting boundary -layer formulation. 
Nevertheless, results obtained with this code reported in Refs. 5 & 8, and in the 
present study indicate that the method yields remarkably accurate results for the 
pressure distribution near shock waves. 

The boundary- layer development on the airfoil and in the wake is determined from 
simple integral methods. The laminar boundary layer starting from the stagnation 
point is computed by an extention of Thwaites 1 method to compressible flow (Ref. 9). 

A transition to turbulent flow is made at a given point on the airfoil on the basis 
of either assigned position or by the use of established transition criteria. The 
method also checks for leading-edge separation and whether it is of the long or short 
bubble type and then assigns transition at the separation point. The development of 
the turbulent boundary layer and wake downstream of transition is determined from the 
lag-entrainment method of Green et aj.. (Ref. 10). 

Calculations carried out with a preliminary version of code (Refs. 5 & 8) 
employed a great deal of numerical smoothing of the surface-pressure distribution and 
boundary- layer parameters near the shock wave and trailing edge. It has since been 
determined (Ref. 11) that the smoothings were responsible for the "spiked" pressure 
distributions observed near the trailing edge in these early studies. We have also 
found that numerical smoothings are not necessary for convergence and they have, 
therefore, been eliminated from the method. 
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There have been a number of other recent attempts to incorporate boundary- layer 
corrections into full-potential flow codes. The work of Bauer, Garabedian and Korn 
(BGK) (Ref. 12) is based on a non-conservative (N-C) treatment of the full potential 
equation and a Nash-McDonald integral method for the turbulent boundary layer. It 
does not solve for the laminar boundary layer near the leading edge and does not use 
transition criteria to start the turbulent flow. It treats boundary- layer displace- 
ment effects on the airfoil through a displacement-thickness formulation, requiring 
repeated conformal mappings of the equivalent airfoil to a circle. Similar methods 
have also been developed by Bavitz (Ref. 13) and Carlson (Ref. 14). These methods 
take no account of wake-induced viscous effects nor strong- interaction effects near 
shock waves or trailing edges. They all employ extensive numerical smoothings which 
are apparently needed to assure convergence of the inviscid/boundary layer iteration. 
These codes also contain adjustable parameters which can be used to improve agreement 
with experiment. The non-conservative formulation employed in these methods leads to 
the generation of spurious mass at the shock wave which causes an extra "sink” drag. 

A post-solution correction for the spurious sink drag was added to a later version of 
the BGK code (Ref. 15) resulting in improved drag predictions. 

A similar method was also developed by Lock and Collyer (Refs. 16, 17 & 18) of 
the Royal Aircraft Establishment (RAE). They employed a more complete boundary- layer 
formulation, similar in many respects to the present method. A surface-source 
formulation was used to represent displacement effects on the airfoil with both 
wake-thickness and wake-curvature terms included in the matching conditions. 
However, strong interaction effects at the trailing edge were not taken into account. 
Numerical difficulties in implementing the wake-purvature condition were reported and 
extensive numerical smoothing of the solution near the trailing edge was required to 
obtain converged solutions. It is likely that the neglect of strong interaction at 
the trailing edge was at least partially responsible for these problems. 

There has been some controversy regarding the proper choice of differencing, 
fully-conservative (F-C) or non-conservative, to be used in these calculations. As 
mentioned previously, non- conservative differencing introduces large errors into the 
drag prediction that require correction. There are also indications (Refs. 5, 8 & 
11) that the pressure rise across the shock wave is underestimated with a non- 
conservative formulation. Lock (Ref. 19) also produced some results with a version 
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of the RAE viscous method employing a quasi-conservation formulation showing shock 
waves that are too far aft. To improve the performance of their method, a partially 
conservative (P-C) scheme was also developed by the RAE (Ref. 18) consisting of a 
linear combination of non-conservative and quasi-conservative differencing. The 
method contained an arbitrary parameter, X, that defines the weighting of the 
non-conservative and quasi-conservative formulae, and which was adjusted to obtain 
good agreement in shock position for a few selected cases. The present study indi- 
cates that the fully-conservative method clearly gives the best predictions of shock 
position, shock strength, and drag provided that the complete boundary- layer formu- 
lation of the present method is employed in the computations. 

It is difficult to form a firm opinion regarding the adequacy of interacting 
boundary- layer theory for the airfoil problem on the basis of results obtained with 
either the BGK or RAE methods. The neglect of potentially important viscous effects 
in the theoretical formulation and the use of extensive numerical smoothing and 
adjustable parameters in the computer code definitely obscure the meaning of the 
results obtained to date. The shortage of reliable experimental data with 
boundary- layer information and the uncertainties associated with wind-tunnel wall 
effects have also impeded progress on the evaluation of the theoretical methods. 
Accurate Navier-Stokes (NS) solutions obtained on fine grids could fill this void 
but, unfortunately, suitable solutions have not yet appeared in the literature. 

The present work is a step in the direction of eliminating some of the above 
problems. For airfoils with cusped trailing edges, the present method includes all 
of the leading-order viscous terms consistent with a rational asymptotic analysis, 
aside from the neglect of normal pressure gradient effects at shock waves. The 
absence of these latter effects in the theoretical model does not appear to be of 
great consequence if the shock wave is not "too close" to the trailing edge. The 
present method also avoids all numerical smoothings, contains no adjustable para- 
meters (aside from those appearing in the turbulence model and transition criteria), 
and yields accurate solutions to the interacting boundary- layer equations 
unencumbered by extraneous numerical issues. The organization of the remainder of 
this paper is outlined below. 
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In Section 2 we describe the viscous-interaction formulation used in the present 
work including a brief description of the classical boundary- layer formulation and an 
outline of the Melnik-Chow strong-interaction trailing-edge theory. In Section 3 we 
describe the boundary value problem governing the outer inviscid flow and outline the 
numerical procedures used to solve the inviscid equations. In Section 4 we discuss 
the methods used to determine the solution in the laminar and turbulent boundary 
layers and to evaluate the matching conditions. In Section 5 results are given for 
various airfoils which illustrate the importance of the individual viscous effects 
that appear in the theoretical model. These include wake thickness, wake curvature, 
and strong interactions at the trailing edge. Comparisons with experimental data are 
also provided including some date for several important boundary- layer quantities. 
In Section 6 we discuss the significance of this work, summarize the principal 
conclusions of the study, and identify related problem areas requiring further 
research. An epilogue is presented in sections which summarizes some improvements 
that have been made to the method since a preliminary version of this report was 
completed in 1980. 
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2, VISCOUS FLOW THEORY 


In this section we describe the viscous flow formulation employed in the present 
work. We will review conventional boundary- layer theory, describe the development of 
the local trailing-edge solution, and will indicate how the conventional matching 
conditions can be corrected for strong interactions at the trailing edge. We first 
briefly review previous theoretical studies of viscous interaction on airfoils. 

2.1 BACKGROUND 

Traditionally, two classes of methods have been used for viscous interactions on 
airfoils at high Reynolds numbers. One is based on the direct numerical solution of 
the full Reynolds equations of turbulent flow (Ref. 20) while the second is based on 
Interacting Boundary Layer Theory (IBLT) (e.g., see Ref. 21), A combination of 
these methods leading to a direct numerical solution of a "thin- layer” or parabolized 
approximation to the full Reynolds equations have also recently appeared (Refs. 22, 
23 & 24). Although methods based on the numerical solution of the full or 
approximate Reynolds equations should be expected to yield the best solutions, they 
have obtained only qualitative results to date. The methods are computationally 
expensive and, consequently, have suffered from poor resolution due to the use of 
coarse meshes - particularly in the streamwise direction. At present, it seems that 
more accurate solutions can be obtained with far less computing time with a 
boundary- layer type approach that takes advantage of the availability of fast 
numerical methods for solving the full-potential equation. These methods permit the 
use of relatively fine grids leading to solutions with good spatial resolution. 
Because of their basic simplicity and high speed, boundary- layer methods have gained 
widespread favor in the engineering community. 

It should be noted, however, that although boundary- layer methods have been used 
for many years and have become one of the standards methods for computing viscous 
flows over airfoils, a completely satisfactory interacting boundary- layer 
formulation for this problem has not yet been achieved. The standard boundary- layer 
approach leads to a description of the flow in terms of coupled inviscid and 
boundary- layer equations which are solved by iteration to obtain self-consistent 
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solutions. The theory is based on an asymptotic expansion of the full viscous 
equations in the limit of Reynolds number. Re «>. The asymptotic theory leads to a 
formal derivation of the matching conditions coupling the inviscid and viscous flows. 
In this way, several terms appear in the matching conditions, all of which should be 
included in a consistent, lowest-order description. Unfortunately, in most previous 
works potentially important terms associated with the wake have been dropped from the 
matching conditions. In addition, there are local failures of the boundary- layer 
approximations in strong-interaction regions near shock waves and trailing edges. 
Thus, a complete interacting boundary- layer formulation should, in principle, 
consider all of the following factors: 

o displacement-thickness effects on the airfoil 
© displacement-thickness effects in the wake 
o wake-curvature effects 

© strong-interaction effects at the trailing edge 
© strong-interaction effects at shock waves. 

The relative magnitude of these effects can be expressed in terms of a basic 

small parameter, s, related to the Reynolds number, Re, such that z vanishes as 

Re ** «>. For turbulent flow it is convenient to identify z with a typical value of 

the friction velocity on the airfoil surface. For the choice z = 0(l/£nRe), the 

boundary- layer thickness, 5, is 0(s) and the displacement and momentum thicknesses, 

2 

6* and 0, are both 0(s ). The boundary- layer formulation can be developed in a 
formal asymptotic expansion in the limit z ** 0. A formal asymptotic theory has been 
developed for non-interacting turbulent boundary levels by Mellor (Ref. 25), Bush and 
Fendell (Ref. 26), Yajnik (Ref. 27), and Afzal (Ref. 28) leading to a boundary layer 
with a law-of -the-wake/law-of -the-wall , two-layer structure. 

The displacement effect arises from a matching of the inviscid and 
boundary- layer solutions. It produces an 0(z ) perturbation to the outer inviscid 
solution which can be incorporated into the solution through the concept of an 
equivalent displacement body or equivalent source flow at the surface. In the first 
approach, the equivalent body is formed by adding the boundary- layer displacement 
thickness to the airfoil and trailing streamline, while in the source formulation, 
sources and sinks are placed on the airfoil and trailing streamline. The equivalence 
of these formulations (and others) was demonstrated by Lighthill (Ref. 29). The 
works of Stuper (Ref. 30) (1933), Pinkerton (Ref. 31) (1936), and Nitzberg (Ref. 32) 
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(1946) -in a rough-way, and Preston (Refs. 33 & 34) (1945, 1949) more accurately, were 
the first to apply the Prandtl boundary- layer theory to compute displacement effects 
on airfoils. Preston used the surface-source formulation to account for displacement 
effects on both the airfoil and wake by distributing sources and sinks along the 
airfoil and wake. These methods have come into common usage although there has been 
a tendency in recent work to ignore, without justification, the displacement effects 
induced by the wake. 

The wake-curvature effect arises from the turning of low momentum fluid in the 

wake along curved streamlines. The turning of the flow induces a pressure drop 

across the curved wake that leads to additional terms in the matching conditions. 

These terms enter the formulation of the outer inviscid problem through a boundary 

condition specifying a discontinuity in pressure along the trailing streamline. The 

2 

magnitude of the effect is 0(s ) and it produces effects on the outer inviscid flow 

that are similar to a jet flap with a negative jet-momentum coefficient of strength 
2 

0(e ). From this is follows that wake curvature leads to a reduction in left that is 
the same order as that produced by the displacement thickness. The importance of the 
wake-curvature effect and its relationship to the jet flap was first noted by Preston 
(Ref. 35) in 1954 (see also Spence and Beasely - Ref. 36). However, serious attempts 
to include wake -curvature effects in interacting boundary- layer calculations were not 
undertaken again until the mid 1970 ! s in a series of papers by Firmin, Hall, Lock, 
and Coyller of the RAE (Refs. 16-19 and 37-39) and by the present authors in Refs. 4, 
5 and 8. Although there were differences in the individual results, all displayed 
significant effects of wake -curvature on section characteristics. 

Interacting boundary- layer theory, as just described, is not uniformly valid at 
trailing edges of airfoils. Non-uniformities in the theory are caused by 
singularities in the inviscid solution for the pressure gradient and streamline 
curvature at the trailing edge and by the change in the no-slip condition across the 
trailing edge. The singularities of the inviscid solution lead to large local 
streamline curvatures and to a breakdown of the boundary- layer approximations near 
trailing edges. Because of this non-uniformity, the pressure can no longer be 
regarded as constant across the boundary layer near trailing edges. This is 
especially important since pressure variations across the boundary layer affect the 
Kutta condition and, therefore, lead to relatively large global effects. Two types 
of singularities arise in the outer inviscid solution. One is related to lift and 
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arises at the trailing edge of any lifting airfoil carrying a non-zero load at the 
trailing edge. The second, due to the stagnation point singularity, appears at the 
trailing edge of an airfoil with a non-zero included angle. 

An analysis of the strong interaction problem associated with lift was initiated 
by the first two authors in Ref. 4 and completed in Ref. 5. In these papers it was 
shown that the inviscid flow leads to square root singularities in the pressure 
gradient and wake curvature at the trailing edge and that, as a result, normal 
pressure gradients were important in the boundary layer near the trailing edge, and 
must be retained in a consistent description of the flow in these regions. It was 
also concluded that a self-consistent boundary- layer formulation and the related 
concept of an equivalent displacement body was not valid near trailing edges in 
turbulent flow. 

A rational asymptotic theory was developed to describe the local flow near 
trailing edges of cusped airfoils. The analysis was carried out for compressible 
flow under the assumption that the flow as subsonic in the immediate vicinity of the 
trailing edge. It was also assumed that the boundary layer was an attached, 
fully-developed turbulent flow near the trailing edge and that the velocity profiles 
could be adequately represented by Cole’s law-of-the-wake/law-of-the-wall form. 
Under these assumptions, it was shown that, in the limit c 0, the boundary layer 
developed a three- layer structure near the trailing edge that is a generalization of 
the two-layer structure arising in non-interacting boundary layers. The three-layer 
structure extended over a streamwise extent of the order of a boundary- layer 
thickness (i.e., 0(e)). The pressure distribution was completely determined from the 
solution in the outermost region, encompassing most of the boundary layer. To 
leading order, the flow in the outer region is governed by linearized inviscid, 
rotational flow equations. The dominant vorticity in the trailing-edge region is 
"frozen 11 along the streamlines of the inviscid flow and is therefore completely 
determined by the upstream velocity profiles. 

The theoretical model accounts for both the pressure variation across the 
boundary layer and for the vorticity in the boundary layer. The two inner layers 
were asymptotically thin compared to the outer, inviscid shear layer and, hence, did 
not effect the latter to lowest order. Thus, the solution of the outer problem for 
the pressure distribution was independent of the Reynolds stresses and closure 
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assumptions except through their influence on the shape of the velocity profile just 

upstream of the trailing-edge region. Analytic solutions were given for the local 

trailing-edge problem in Ref. 5. These solutions showed that the pressure 

3/2 

perturbations generated in the trailing edge region are 0(e ) which is an order of 

2 

magnitude larger than the 0(s ) disturbances generated by conventional displacement 
and wake-curvature effects. 

Similar interaction effects are also expected from the non-uniformity associated 
with the stagnation point singularity arising on airfoils with non-zero trailing-edge 
angles. Theoretical evidence indicates that the turbulent interaction changes the 
pressure by 0(1) from the stagnation pressure arising in the inviscid solution. 
Thus, this effect is clearly larger than the others considered so far and can be 
expected to have important effects on the drag of airfoils with wedge-shaped trailing 
edges. Unfortunately, definitive studies of this important strong-interaction 
problem have not yet been carried out. 

Aside from the non-uniformities stemming from singularities in the inviscid 
pressure distribution, an additional non-uniformity arises from the discontinuous 
change of the no-slip boundary condition across the trailing edge. The jump in 
boundary condition at the trailing edge causes an additional breakdown of the 
interacting boundary- layer solution. A similar effect arises in laminar flow with 
the transition of the laminar boundary layer to a wake across the trailing edge. In 
this case, local asymptotic solutions were developed for the non-lifting flat plate 
by Stewartson (Ref. 40), Jobe and Burggraf (Ref. 41), Veldman and van de Vooren (Ref. 
42), and Chow and Melnik (Ref. 44) and for the lifting case by Brown and Stewartson 
(Ref. 43) and Chow and Melnik (Ref. 44). These solutions were based on the 
triple-deck theories of laminar interacting flows. Unfortunately, because of the 
difference in structure of laminar and turbulent flows, the laminar triple-deck 
theories do not apply to the turbulent problem. In turbulent flow the transition of 
the boundary layer to a wake is largely controlled by the law-of-the-wall region 
which has no counterpart in laminar flow. Since the law-of-the-wall region is 
exponentially thin compared to the boundary- layer thickness, the transition process 
will have no effect on the leading-order outer solution for the pressure distribution 
in turbulent flow, while in laminar flow it has a dominant effect. In general, the 
turbulent boundary layer/wake transition process in the wall layer will have only a 
higher-order effect on the outer solution and need not be considered in the solution 
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for the pressure distribution and section characteristics. However, for a symmetric, 
cusped airfoil at zero incidence the pressure-induced, strong interaction vanishes 
and the local flow is dominated by the turbulent boundary/ layer/wake transition which 
then becomes the most important interaction effect near the trailing edge. Although 
there have been some exploratory studies of this latter problem by Robinson (Ref. 45) 
in 1967 and more recently by Alber (Ref. 46), these works were not definitive and 
they did not establish the correct asymptotic flow structure at large Reynolds 
numbers. Fortunately, the lack of knowledge of this aspect of the solution will only 
affect higher-order terms and need not be taken into account in order to determine 
the leading-order solution for the pressure, lift and drag. 

For transonic flow an additional strong interaction arises at points where shock 
waves impinge on the boundary layer. In turbulent flow, discontinuities in pressure 
across the shock wave induce corresponding discontinuities in displacement thickness 
(Refs. 47-48), leading to a breakdown of theories based on an interacting 
boundary- layer description. The shock wave penetrates into the boundary layer and 
generates large pressure variations across the layer that invalidate the standard 
boundary layer approximations. 

Rational asymptotic theories, recently developed (Refs. 47-53) to describe 
shock-wave/turbulent boundary- layer interactions, are based on a large Reynolds 
number (e 0) asymptotic expansion of the full Reynolds equations. The boundary 
layer develops a multi-layer structure, similar to the turbulent trailing-edge 
interaction discussed previously. The pressure disturbances generated by the 
interaction were shown to be the order of e x (shock strength) and hence larger than 
the classical displacement effect for shock strengths greater than 0 (e). The 
analyses in Refs. 47-53 considered weak shock waves of various strengths. It was 
also demonstrated in these works that the shock strength required to separate a 
turbulent boundary layer is 0(1). These works considered only local behavior and did 
not address the question of integrating these solutions into a complete airfoil 
problem. Fortunately, this is not a major concern here since the normal pressure 
gradient effect in shock-wave/boundary- layer interactions seem to have only a highly 
localized effect on the pressure distribution. Alternative formulations of 
shock-wave/boundary- layer interactions in the turbulent flow have also been proposed 
recently by Inger and his co-workers (Refs. 54 & 55). This work was based on 
application of Lighthill*s laminar interaction theory (Ref. 56) to the turbulent 
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problem. Inger's theory represents an ad-hoc approach to the problem that is not 
related to a rational asymptotic description of the flow. 

2.2 INTERACTING BOUNDARY -LAYER THEORY 

We consider the problem of a uniform viscous flow over an airfoil of chord, c, 
free-stream speed, U , and density, p , at an angle of attack, a, and Mach number, 

oo oo 

M . The main assumptions of the work are that the Reynolds number Re is large in an 
asymptotic sense and the free-stream Mach number is less than one. Under these 
conditions the problem can be analyzed by rational asymptotic methods based on 
expansions of the full Reynolds equations in the limit Re 4 or 8 ^ 0. In this 
limit, a boundary- layer type flow structure develops as sketched in Fig. 1. The flow 
outside the strong- interaction region can be described by a standard 
inviscid/boundary- layer formulation with conventional matching conditions coupling 
the solution. In the standard approach, the solution in the outer region is 
represented by a sum-type asymptotic expansion. The solution in the outer region is 
inviscid to all orders, apart from exponentially small terms, with the leading term 
given by a solution to the full non-linear inviscid equations. Second and 
higher-order terms of the outer solution are governed by linearized inviscid flow 
equations . 


This representation is somewhat inconvenient for supercritical flows because of 
the need to deal with perturbations of discontinuous solutions. This problem can be 
avoided by use of a slightly different representation for the outer solution. In 


EXTERNAL INVISCID FLOW 



Figure 1 Flow Field Regions at High Reynolds Numbers 
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this alternative formulation, the outer inviscid flow is governed to all orders, by 
the full non-linear inviscid-f low equations with boundary conditions on the airfoil 
and in the wake determined from the viscous matching conditions. Higher-order 
viscous effects on the outer inviscid flow appear only through the matching 
conditions. Thus, in this formulation the full non-linear inviscid equations must be 
solved repetitively, employing successively improved matching conditions. 

The form of the solution in the inner boundary- layer region depends on whether 
the boundary layer is a laminar or turbulent. In the laminar case, the solution has 
the form given by higher-order boundary- layer theory with the leading term governed 
by the usual non-linear Prandtl boundary- layer equations. For turbulent boundary 
layers, the form of the solution changes because of the two-layer structure of 
turbulent boundary layers at high Reynolds numbers. The appropriate asymptotic 
theory has been worked out for fully-developed turbulent boundary layers in Refs. 
25-28. The leading-order solution is governed by boundary- layer equations with 
linearized convective terms. However, the usual non-linear boundary- layer equations 
do provide a composite equation that is valid to second-order in both the wake and 
wall regions of the turbulent boundary layer. In the present work we follow the 
composite equation approach and use the standard non-linear boundary- layer equations 
to describe the flow in both the laminar and turbulent boundary- layer regions. 

The outer inviscid and inner boundary- layer solutions can be combined in the 
usual fashion to form a single composite expression that is uniformly valid in both 
regions. We employ curvilinear coordinates (s,n) and corresponding velocity 
components (Tl,!)) along and normal to the streamline defining the airfoil and wake. 
In this study we employ nondimensionalized quantities with the velocities and density 
normalized by free stream velocity and density (U , p ) , the pressure and shear 

OO CO 

2 

stresses by twice the dynamic pressure, (p U ) and all lengths are scaled by the 

OO OO 

chord (c) of the airfoil. Composite expressions are denoted by the upper case script 
letters. The composite solution for the velocity components and the pressure are 
then represented in the form of a sum of the inviscid solution (capital letters) plus 
boundary layer solution (lower-case letters) minus the common part (subscript cp) as 
follows : 


'll =U(s, n; v 0 ; JV]]; r) + u(s, n) -u cp (s, n) 


(la) 
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U = V (s, n; v 0 ; ([VD; D+v(s, n)-v cp (s, n) 


(lb) 


(P = P(s, n; v 0 ; [[VD ; D + p(s, n)-p cp (s, n) 


(lc) 


where n is a stretched inner variable defined by 

n = n/e (2) 

Similar expression are also assumed for the density, (R , and temperature, T • The 
quantities, v^, j£vj] and T appearing in the inviscid solution are functions entering 
the boundary conditions that are determined by matching the inviscid and 
boundary- layer solutions. The function v^(s) specifies the source velocity on the 
surface of the airfoil, [£v(s)J] specifies the jump in source velocity across the 
wake and T(s) is the circulation distribution along the wake. For an inviscid flow 
all three functions are zero as are the contributions from the boundary- layer and 
common-part term in the above representation. The solution is then given by the 
first term with boundary condition V(s,o) = 0. 

If u(s,n) and p(s,n) are the boundary- layer solutions for the streamwise 
velocity and density profiles, the boundary- layer solution for the normal component 
of velocity and pressure distribution can be written in 

v(s, H)=v(s ) 0) - dPe . Ufl n+ — -j— f (p e U e -pu)dn + 0(e 3 ) (3a) 

P e as p e ds Jo 

p (s, n) = p (s, 0) — /c(s) p e U p n+jc(s) f (p e U 3 -pu 2 ) dn + 0 (£ 3 ) (3b) 


where k(s) is curvature of the airfoil or the wake streamline and where and are 
the surface values of the inviscid solution defined by 

p e (s) = limR (s, n) (4a) 

n- 0 

U e (s) = lim U(s, n) (4b) 

n- 0 

We introduce here the definitions of the surface values of the normal velocity and 
pressure from the outer inviscid solution as: 

V e (s) = limV(s, n) (4c) 

n - 0 
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P e (s) = lim P(s, n) 

n - 0 

(4d) 

Note that the above quantities are generally different for n 0 + , the upper surface 
and n 0 , the lower surface. When necessary to distinguish between the surfaces we 
will include a + or - as a superscript. 

Comparisons of the boundary- layer solution for n -> <» to the 
for n 0 leads to the following identification of the common parts: 

inviscid solution 

Ujp — U e (s) +0(e 2 ) 

(5a) 

v cp = v(s, 0) - e n + e 2 v„ (s) +0 (e 3 ) 

P e as 

(5b) 

p cp = p(s, 0)- k(s) p e Ug n + e 2 p CT (s) +0(e 3 ) 

(5c) 

where 


, , . d p p Up6* 

” (s)= 

(6a) 

p„(s)=p e Ue «(s) (<5*+0) 

(6b) 

and where 5* and 0 are the usual boundary- layer displacement and momentum 
thicknesses. Note that 6* is defined such that it is positive on both upper and 
lower surfaces. Matching conditions on the airfoil surface can.be derived by 
substituting the above expressions for the common parts into Eq. (1). The 

requirement that the normal component of velocity , *U , vanish on the airfoil surface 
then leads to the usual matching condition on V(s,n) for n 0 : 

V e (s) = e 2 v 0 (s) + 0(e 3 ) 

(7) 

Similarly, substitution of Eq. (5c) into Eq. (lc) leads to the following expression 
for the surface pressure distribution 

(P e (s) = (P (s, o) =P e (s) -e 2 p 0 (s) +0(e 3 ) 

(8) 


where P^Cs) is the surface value of pressure as determined from the inviscid solution 
and the second term in Eq. (8) represents a correction for the pressure variation 
across the boundary layer. 
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The corresponding expressions in the wake are somewhat more complicated because 
of discontinuities of the common-part terms across the wake. If we impose the 
condition that the full solution forTljTJjCP, etc. is continuous across the wake and 
take proper care for the discontinuities of the common-part terms we can arrive at 
the following jump conditions for the outer inviscid solution: 

IVD =V(s, 0 + ) — V(s, O')=e 2 v w (s)+O(e 3 ) (9a) 

EP]]=P(s, 0*) — P(s, 0-) = -e 2 p w U w c w K(s)+O(€ 3 ) (9b) 

Moreover, since 

Euu = -r-^r EpB +o([[p] 2 ) 

then 

ffUB =U(s, 0 + ) — U(s, 0")=e 2 c w K(s)+O(€ 3 ) (9c) 

where 

e 2 v w (s) = v;(s) - v- a (s) = - d - ^ ^ (9d ) 

e 2 c w (s) = — U w (<5 * + 0 w ) (9e) 

and where, p and U are the surface values of the outer inviscid solution on the 
w w 

JU 

wake streamline, and 6^ and 0^ are the respective sums of the upper and lower 
displacement and momentum thicknesses. The composite solution for the pressure at 
the wake axis is continuous and is given by 

(p e (s) = (p (s, o) = p+ (s) - p;(s) = p;(s) - P ;(s) (10) 

where the plus and minus superscripts refer to the upper and lower sides of the wake 
axis. (In deriving the above expressions we have used the fact that the jumps in 

2 

and across the wake axis are 0(e )). Equations (9a) and (9b) are jump conditions 

to be imposed along the wake axis of the outer inviscid flow. It should be noted 

that Eqs. (9) imply that the streamline slopes of the outer inviscid flow are also 

2 

discontinuous across the wake axis to 0(e ) and that this in turn implies that the 
wake axis cannot be chosen as a streamline of the outer inviscid solution (to this 
order) . A convenient choice for the wake axis is the composite solution for the 
streamline passing through the trailing edge. An equivalent analysis of the matching 
conditions in the wake is included in the more general work of Viviand (Ref. 57). In 
the present study we use an irrotational approximation to justify the introduction of 
a velocity potential in the outer inviscid problem. The above matching condition on 
the streamwise velocity component implies the velocity potential, must also be 
discontinuous across the wake. The appropriate jump condition on $ can be determined 
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from Eq. (9c), rewritten in the form 


dr 2 / \ 

— =e'c w K(s) dla) 

where the circulation strength, T, is equal to the jump in potential across the axis. 
That is , 

r= $(s, o + ) -$(s, o") (iib) 

If B is defined as the angle between the wake streamline (defined by the composite 
solution) and the airfoil chord then, 


= d£ 

ds 


( 12 ) 


The expression for T can then be written in differential form, 

dr=e 2 c w d/3 (13) 

relating changes in circulation to the turning angle of the wake. The interested 
reader is referred to Refs. 5, 8, 16, 17, and 39 for further discussion of this 
condition . 


The complete set of matching conditions derived above are summarized in Fig. 2. 
Note that within this formulation the wake thickness is prescribed but the location 
is not. The location of the wake is determined as part of the outer inviscid 
solution and is free to assume an equilibrium position consistent with the prescribed 
pressure jump across the wake. 



Figure 2 Viscous Matching Conditions 

For cusped airfoils at lift, both the pressure gradient and streamline curvature 
generally become unbounded at the trailing edge with the following behavior: 
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(14b) 


lim 


8 - 8 


te 



s-s te l 


-i/2 


where the subscripts (te) denote values at the trailing edge, R and U^ e are the 
density and speed at the trailing edge, is a constant that depends on the overall 

inviscid solution, B is the Prandtl-Glauert factor - M^, and the upper and lower 
signs in Eq. (14a) refer to the two sides of the airfoil. From the momentum integral 
equation and the definition of the surface-source velocity in Eq. (6a), we can 

show that Eq. (14a) implies a corresponding singularity in v^, as follows 


lim 

8 - 8 * 


V a =± 


11 * 0*61 


te 


(^r) ls - ! 


*te 1 


- 1/2 


(Uc) 


where 6 + is the boundary- layer displacement thickness of either the upper or lower 
surface boundary layers at the trailing edge. In deriving this result we assumed the 
total temperature is constant across the boundary layer and have used the large 
Reynolds number limit for the shape factor, h. That is, 


h= h te = 1 + (r — 1) Mj e for R e — « 


Similarly, from Eqs. (9b), (9e) and (14b) we obtain: 

lim g PB = Pfft U * C< x (< y I s - s te I ~ 1/2 (14d) 

®" s te 

Detailed analysis indicates that these singularities lead to a breakdown of 
conventional boundary- layer theory at trailing edges and to a growth in the magnitude 
of viscous effects induced by both the displacement thickness and wake curvature 
terms from O(e^) to 0(e^^) near trailing edges. A separate "inner” solution is 
required to resolve the nonuniformity at the trailing edge as described in the next 
subsection. 


2.3 TRAILING-EDGE REGION 

In this section we present the results of an analysis of the strong-interaction 
region that develops in turbulent flows near trailing edges. An analysis for the 
incompressible case was first presented in Ref. 4 and later extended to compressible 
flow in Ref. 5. The analysis is based on a formal application of the method of 
matched asymptotic expansions (Ref. 58) to the full Reynolds equations of turbulent 
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flow for Reynolds number Re ** «>. The theory was developed under the following 
assumptions : 

o the trailing edge of the airfoil is cusped 
o the flow is locally subsonic in the trailing-edge region 
e total temperature is constant in the trailing-edge region 

e the boundary layer upstream of the trailing-edge region is a fully-developed 
turbulent flow 

o the boundary layer approaching the trailing-edge region is not separated 
In the present context, the assumption of a fully-developed turbulent flow is taken 
to mean that the velocity profile in the outer part of the boundary layer can be 
adequately represented by a Coles* law-of-the-wall/ law-of-the-wake, with a small 
velocity defect. That is, we assume the velocity profile in the outer part of the 
boundary layer, upstream of the trailing edge can be represented in the form 

u = u e (s)[l+€f(s, n) + **-] 

with 
and 

U* = \/t w /p 6 (15c) 

where x^, u*, and 6 are the local skin friction, friction velocity and boundary- layer 

thickness, k is the von Karman constant, it is Coles* wake parameter, and W is the 

wake function which, in the present work, is assumed to be expressible as a simple 

third-order polynomial. It is also assumed that the profile function, f(s, n) 

approaches a definite limit, f (n) at the trailing edge, s = s . The above 

re t e 

assumptions permit a relatively complete solution to be obtained for the local flow 
near the trailing edge. In this region the solution develops the multi-layer 
structure shown schematically in Fig. 3. The upstream flow is divided into 
conventional inviscid and boundary layer regions over a streamwise length of 0(1). 
The turbulent boundary layer has a two-layer structure consisting of an outer, wake 
region and an inner, equilibrium wall layer. The velocity profile in the outer 
region is described by Eq. (15). 


(15a) 

(15b) 


The velocity profile in the inner layer is expressed 
which for incompressible flow is written as 


u = eU e (s)F 



in a law-of-the-wall form. 


(16) 
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€ = V C fo /2 (FRICTION VELOCITY) ~ 0 (LOG Re) 
? = {e 2 Re)’ 1 - 0{L0G Re)/Re 


Figure 3 Flow Field Structure Near the Trailing Edge 

where s is the ratio of wall-layer to boundary- layer thickness and is given by 

e=(e 2 R e )-1 (17) 

It follows that the wall layer is exponentially thin compared to the boundary layer 
thickness . 

The formal justification of the two- layer structure in terms of an asymptotic 
limit solution for Re ** <*> was carried out by Mellow (Ref. 25), Bush and Fendall (Ref. 
26), and Yajnik (Ref. 27) for incompressible flow. These studies showed that the 
velocity profile functions satisfied simplified boundary layer equations with 
linearized convective - terms . The usual non-linear boundary layer equations can be 
viewed as composite equations valid in both inner and outer regions and containing 
additional higher order terms. 

An extension to compressible flow was attempted by Afzal (Ref. 28) but the 
analysis was in error, with the result that his assumed inner and outer expansions 
did not match. The reasons for the error and the modifications necessary to obtain a 
valid compressible solution are discussed in Ref. 47. It was shown there that the 
small defect form given in Eq. (15) carries over to the compressible case, but that 
the inner expansion required modification to incorporate a preliminary van Driest 
compressibility transformation. One of the more firm rational results in turbulence 
theory concerns the logarithmic behavior of the velocity profiles in the overlap 
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region between the inner and outer regions, i.e., n/e << 1 and n/ee >> 1 leading to 
the form of the velocity profile given in Eq. (15). 

The present work can be viewed as an extension of the asymptotic theory of 
non-interacting boundary layers to the strong interaction problem. The boundary 
layer in the interaction region develops a compact, multi-layer structure as 
sketched in Fig. 3. Three layers are required: 

1. An outer layer that is inviscid and rotational. The vorticity arises 
from the nonuniform flow in the upstream boundary layer. In a fully-developed 
turbulent flow the total vorticity is small and the solution in this region 
can be described as a perturbation of the external potential flow induced by 
the small vorticity. To lowest order, it can be assumed that the vorticity 

is convected along the nearly parallel streamlines of the potential flow. In 
the terminology of secondary-flow analysis, the flow in this region is class- 
ified as a small shear, small-disturbance flow (Ref. 59). 

2. An inner layer next to the wall. In this region the flow is de- 
scribed as equilibrium wall layer that is a continuation of the wall layer 
from the upstream flow. The total stress, laminar plus Reynolds stress, is 
constant across the layer. The thickness of the layer is exponentially small 
compared to the outer layer. 

3. A blending layer situated between the outer and wall layers. It 
is necessary in order to match the Reynolds stresses in the outer and wall 
layers. The layer is thinner than the overall boundary layer but thicker 
than the wall layer. 

Turbulent closure approximations are required to -lowest order, only in 
the two inner layers. The inner layers are thin and do not affect the 
determination of the solution in the outer inviscid region. The outer 
region is most important since it alone determines the pressure distribution 

Outer Exp ansi on and Governing Equations . The solution in the outer region is 

developed as a perturbation to the basic external inviscid flow in the following 
form: 

u = U(x, y;e) +e U te f(y)+e 3/2 u (1 >( x , y) + * • • (18) 

v=V (x, y;e)+e 3/2 v (1) (x, y)+*** ( 19 ) 
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p=P(x, y;e) +e 3/2 p ( 1 ) (x, y)+*** 


( 20 ) 


where (x,y) form a local Cartesian coordinate system with an origin at the trailing 


edge and scaled as follows 
X = S g ~^ Ste » y = n/€ 

with the Prandtl-Glauert factor B defined as 


B= v/l-Mf e 


( 21 ) 

( 22 ) 


The leading terms of Eqs. (18) - (20) are meant to include all terms that result 
from expanding the known inviscid solution in powers of (s,n), transforming to (x,y) 
and expressing the result as an expansion in z . For a cusped airfoil this leads to 
the expansion 

U = U te [l+e 1 / 2 U 1 (x, y)/B +• * • ] (23a) 

V = U te [e 1 /z V 1 (x, y)+*“] (23b) 

where U. , V. are known algebraic functions and U is the trailing-edge value of the 

inviscid velocity. The second term in Eq. (18), f(y), arises from the nonuniform 

flow in the upstream boundary layer and is a known function of the form given in Eq. 
(15). 


Notice that the expansion for u is similar to the law-of -the-wake in the 
upstream boundary layer, except that now the full y-dependence of the inviscid flow 
is accounted for in the leading term, U(x,y;c). The equations governing the 
perturbation terms u^^, p^^ are arrived at by substituting the above 

expansion into the full Reynolds equations of turbulent flow. The series solution 
must be augmented with similar expansions for the density and temperature. For 
convenience we assume the total temperature in the upstream boundary layer is 
constant across the boundary layer. The density can then be computed from the 
equation of state and the total temperature. This is known to be a highly 
satisfactory approximation in the speed range of interest. 


The above analysis leads to a set of partial differential equations governing 
the disturbance to a weakly sheared, compressible, inviscid flow. The disturbance 
equations can be reduced to incompressible form by a generalization of the 
Prandtl-Glauert transformation of subsonic potential flow. After transforming to the 
scaled coordinates x,y by Eq. (21) we transform the dependent variables u , v , 
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u (1 } = U te B” 1 ' 2 [a 2 u + (a 0 - yM \ e ) ^ ^ - ao f U t + ajX 
v <1 > = U te B 1/2 |a 2 v + (l-a 2 )fV 1 -a 1 x 
p ( 1 ) =p te U 2 te B- 1 / 2 ^a 2 p-a 1 x — 

where U^, and are normalized velocities and streamfunction determined from the 
inviscid expansion given in Eqs. (23). For cusped airfoils we have 

Ui — iVj = -i c a (x + iy) 1 / 2 (24d) 

or in real form, 

U 1 = c a r 1/2 sin^ , V A = c a r 1/2 cos-| 


(24a) 

(24b) 

(24c) 


2 3/9 3ft 

*l = -3 c a r 3/2 cos — 

where is a known constant and (r,52) are polar coordinates with r 2 = 
tan J2 = y/x. The function x is the streamfunction corresponding to f(y), 


(24e) 

x 2 + y 2 and 


: = f %) dy 
Jo 


(25) 


The parameters a^, a^, a ^ and are Mach number dependent constants defined by 


a 0 = i(Co-l)-^(Co + l> 

2 

ai= -iP MCo+1) 


(26a) 

(26b) 


a 2 =i(l + B' 2 )(C 0 + l) 

c 0 =i+(r-DM 2 te 


(26c) 

(26d) 


With these transformations the basic equations governing the disturbances in the 
outer region reduce to the following linearized, incompressible equations of mass and 
vorticity 


9£.9v =0 (27) 

9x 9y 


9u 

3y 


9v 

9x 


=*i(x, y) 


d 2 f(y) 
dy 2 


(28) 
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The transformed static pressure is determined from the linearized Bernoulli 
equation 


_ _ rrr df 

p=-u-fu 1 + — 

dy 1 


(29) 


If the disturbance streamfunction, ip, is introduced, Eq. (28) becomes 

o— r]‘- f 

v *=-t = *^(x,y> ( 30 ) 

where u = 8i/;/8 y, v = - 8i^/3x. This is the same basic equation derived in Ref. 4 for 

the incompressible trailing edge problem. Equation (30) is a simple Poisson equation 
relating the disturbance streamfunction, ip, to the perturbation vorticity, £. The 
vorticity perturbation arises from generation in the upstream boundary layer and 
convection along the curved streamlines of the external potential flow, as sketched 
in Fig. 4. The vorticity is known in terms of the upstream boundary layer velocity 
profile and the inviscid streamfunction. The Poisson equation must be solved subject 
to the boundary condition that the normal component of velocity vanish on the airfoil 
surface, which, to the same order of approximation, is represented by a slit along 
the negative x axis. Thus, 


v(x, 0) = 0 for x<0 


(31) 



YM = $*o (Y + e % * 1 + •••) + •• • 

= r 0 (v) + e 1 / 2 r 0 (Y )^ 1 + . . . 

= -ef'(Y)-e 3/2 f"(Y)* 1 + ... 

INITIAL 

PROFILE INTERACTION 

Figure 4 Vorticity Distribution in the Trailing Edge Region 
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In addition, we must impose the condition that the pressure and flow deflection be 
continuous across interfaces at the boundary layer edges and wake centerline. We 
must also assure that the disturbances decay in the far field rapidly enough to allow 
matching to the outer inviscid potential solution. 


It is convenient to introduce the pressure and a variable w, related to the 
disturbance flow angle, as dependent variables. Thus, we transform from (u, v) to 
(p, w) defined by 

v=£ + fVj (32a) 


u= "[ p+fUr | 4, ‘ 


(32b) 


and arrive at the symmetric pair of differential equations 


9p _ jto _ _ 2 r jEi 

9x 9y 9x 


9 p du) 
9y + 9 x 



(33a) 

(33b) 


The function f(y) contains contributions from boundary layers on the upper and lower 
surfaces of the airfoil. Since the equations are linear the two layers can be 
treated separately and the results superposed to obtain a complete solution. 


Scale Transformation . - The solution to the outer problem can be reduced to the 
determination of universal functions through the introduction of the following scale 
transformations appropriate to either the upper or lower surface boundary layer. 
First, change coordinates according to 

77 = y/<5 £=x/<5 v = v/6 (34) 


where 6 is the thickness of the boundary layer at the trailing edge. Then the 
boundary- layer profile functions f and X are expressed in the law-of-the-wake form as 
f = y* [ti t (ri) + TT h' 2 (r))] 0 < 77 < 1 (35a) 

X =7*6 [h^rj) + 7r h 2 (r]) ] 0 < 7j < 1 ^35b) 


h^Tjlogrj -7), h 2 =- f W ( 77 ) dJ 7 

Jo 


(35c) 


and are zero outside the strip 0 < q < 1 (primes denote 3/3 ti). tt is Coles' wake 
parameter, defined by the initial boundary- layer profile, and Y* is a friction 
parameter defined by 
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( 36 ) 


-v* - V^w/ Pte^te 

y - k 

and is the skin friction of the upper surface boundary layer at the trailing edge. 
The following additional scalings are employed 


U t =:c a 6 1/2 U(4, 7J), V 1 = c a 6 1/2 V(^, n) (37a) 

where 

U = P 1/2 sin-^- , V= y 1/2 cos — (37b) 

Li Z 

p» 2 = | 2 +r] 2 , f2 = tan _1 (77/^ ) (37c) 

and 

p = c a Y*«5 1/2 [Ptd, T]) + 7T p 2 (|, T])] ( 38a ) 

u = c a y*6 U2 [ajjd, rj) + 7rw 2 (5» tj)] (38b) 


The functions p^ ^ an< ^ ^ are universal functions of (£, ti) that are independent of 
all parameters appearing in the problem; although p^ and depend on the functional 
form used to represent the wake function W(t}). 


A A 

The p and to functions satisfy the following inhomogeneous Cauchy -Riemann 
equations , 


^-^■= 2h;,jto> BUl £ ^ 

(39a) 

if^ + a|^_ 2hbw %* 

(39b) 

with 


" 1 , 2(5 1 0) =0 for | <0 

(39c) 

Particular Solution. - A particular solution to the inhomogeneous Cauchy -Riemann 
equations can be found by standard complex variable analysis. Following Ref. 60 we 

introduce the complex representation 


z = £ + it] 

(40a) 

U -iV = -iz 1/2 

(40b) 

A l,2 = Pl,2 + i "l,2 

(40c) 


With (z,ti) considered as independent complex variables, Eq. (30) can be written as 
the single complex equation 
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(41) 


0A. < o hj o(7)) 

dr) ~ Vz - 2 \r) 

where the derivative with respect to ti is applied holding z constant. The solution 
to Eq. (41) can be expressed as the sum, 

^1,2 "^pi, 2( z > 7 )) + ^H1,2( z ) ^ 2 ) 

where the first term is any particular solution satisfying Eq. (41) and the second is 

an arbitrary analytic function that is determined from boundary and interface 

conditions. The interface conditions now simply require continuity of X- across 

i , z 

interfaces. A particular integral is easily found by quadrature 


V.2 aFl .2 (Z ’ n)= < 



Vz-2i7} 


7 ] > 1 
0 < 7] — 1 
7) < 0 


The above solution is continuous across the interface, 
discontinuities across the real axis, ti = 0. 


T1 


(43) 

1, but has jump 


The component corresponding to the logarithmic term in the velocity profile is 
found by direct integration. The quadrature for the second component can be carried 
out for polynomial approximations to the wake function W(ti). We employ the 
representation 




(44) 


With Cj = 15.369157 and = -0.36536. Eq. (44) is an accurate (to four places), 

rational approximation to Cole's wake function, 1 + cos itn . Thus, Eq. (44) should 

provide an acceptable fit to most turbulent velocity profiles. These procedures lead 

to explicit algebraic expressions for the particular solution X (z,n). They are 

px , z 

somewhat lengthy and, hence, will not be written-out here. For the purposes of the 

present study we require only the limiting values on the real axis, ti -*• 0 + , leading 

to the jump discontinuity referred to above. We introduce the real functions 8 (£1 

1,2 

and defined as the real and imaginary values of X ^ ^ as the real axis is 

approached from above. That is, 

0l, 2 (!) +i(T l,2(£) = limX Pl,2( z » V) 
rt- 0 + 

These functions can be conveniently computed from the limiting values of analytic 

functions g (z) determined from X , Thus, 

^ Dl . 2 . 


/ ? l,2 +io 'l,2= lim gl,2( z ) 

»l-0_ 


(45) 
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where 


gj(z) = 2i |vz - V z - 2i + 77 log ^ ' — + 2 -/ z — j ( 46 ) 

g 2 (z) = - ^3- | 1 _ ^-(z -i) 2 j j^[z -i] [z 5/2 -(z -2i) 5/2 ]--~ [z 7/2 - (z - 2i) 7/2 ] 

[ z-i ] U 5/2 - ( z- 2 0 5/2 ] -i[z 5/2 + (z-2i) 5/2 ]j-2iz 1/2 (/ 

The only singularities in the strip are branch points at infinity and the origin. 
The appropriate branch of the multivalued functions are defined by a cut along the 
negative real axis and the limit n 0 . The functions 0 (£) and P (£) can be 

1 , Z 1 , Z 

shown to possess the following symmetry 


ff l,2«)=0i,2(-*) (48) 

Thus, only values on the negative (or positive) real axis need be computed. The 
functions are plotted in Figs. 5 and 6. For £ - «> they possess the following 
asymptotic behavior 


^.2 — Ul- I/2 

(49a) 

ft— 

(49b) 

The behavior for £ 0 is given by 


Ol 2 + f +0(4) 

(49c) 

a 2 1.3740 + 0(4) 

(49d) 

/ 3 1 -- 2 -/ry ini 241 +2Vrj+ 0 (4) 

(49e) 

p 2 1.3740 + 2/^f + 0(4) 

(49f) 



Figure 5 Particular Solutions for the Downwash Functions cr-j, 02 
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Figure 6 Particular Solutions for the Loading Functions 

The particular integral specified above is not unique and to it we can add 

arbitrary analytic functions of z to obtain other solutions. From the above 

expressions we can determine the contributions of the particular solution to the 

normal component of velocity on the airfoil and to the jump in pressure across the 

wake. These quantities have a direct correspondence to the source velocity and 

pressure jump appearing in the matching conditions of the external flow discussed in 

Section 2.2. It is convenient to arrange the particular solution so that these 

quantities match the limiting behavior of the external solution for x 0, z 0 as 

given in Eas. (7 and 9). This can be accomplished by adding the following piecewise 

analytic function to the particular solution F (z, n), 

1,2 


AF U2( z ) = (z ~fjT77 V > 0 


where a is a constant. The function AF (z) is analytic in the upper half plane 

o 1,2 

and is defined to be zero in the lower half plane. Hence, it contributes to the 
discontinuity of the particular solution across the real axis. The contribution to 
the jump across n = 0 is given by 


A/? li2 + iA<r lf2 = 


(F+i v 


- jrfa- for U I — 


The full discontinuity of the particular solution across the real axis is therefore 
equal to 
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(52) 


£l,2(£) +icr l f 2(£) ~^l f 2 + icr l f 2 + 



The contributions of the above solution to the normal component of velocity 
on the airfoil and to the pressure jump in the wake are given by 


where 


3 / 2 ; ^ /2 Ar. 

p (1 + tt; 1_ 

7 '(e6B U ) + luJ *( 

S S t e\”| 
e6B )\ 

3/2 n e 3/2 A (" 

p ' =± d^) L 



3/2 A = c a U te B 1/2 

a. 2 y* 6 1/2 (i + 7 r) 



(53a) 

(53b) 

(53c) 


and where v^ and p^ are contributions of the particular integral to the inner 
solution Eqs. (24). The above expressions correspond to both upper and lower 
surfaces with the appropriate values of 6 , Y*, and it. The plus and minus signs in 
Eq. (53b) refer to the upper and lower surfaces, respectively. 


To demonstrate matching to the external solution, we expand the above 
expressions for | £ | ~ using Eqs. (49) and (51) to obtain 


lim .3/2 - _ e 3/2 A(l + a3 ) _ e 2 A(l+ a ,)B 1/2 <5 1/2 

!i” e v " hF^ (s u -s)<« — 

lim _ £ 3/; A(l-aj 

P tt) Ul (s-S te )l/2 


(54a) 

(54b) 


For constant total temperature and large Reynolds number the displacement thickness 
is given by 

5 *-*C 0 (l + 7f)(T*6) for R e — *• °° (55) 

From the above relations and the definition of a^ in Eq. (26c) we can show that Eqs. 
(54) reduce to the singular expressions given in Eqs. (14c, d) provided the constant 
a^ is given the value 

a - — M te 

3 2 •— M(e (56) 

Note that is is the negative of the particular integral that corresponds to the 
surface source velocity in Eq. (7). These results establish the close 
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correspondences between the matching conditions of the external solution and the 
discontinuities of the particular solution on the real axis. The determination of 
the particular integral is the most important result of this section. 

The particular integral provides the contributions generated by the 
inhomogeneous forcing term of the differential equations. It is not a complete 
solution to the problem because it is discontinuous across the wake axis and it has a 
nonzero imaginary part on the negative real axis and, therefore, violates the 
boundary conditions given by Eq. (39c). 

The complete solution in the trailing-edge region can be obtained by adding a 

complementary solution to the particular integral determined above. The 

complementary solution can be represented by an arbitrary analytic function of the 

complex variable, z, which can be adjusted to satisfy the boundary conditions. The 

analytic function is determined from conditions that v = v + v is zero on the 

P “ 

airfoil and the solution is continuous across the wake axis. 

The solutions just described account for contributions generated by the upper 
surface boundary layers. Similar contribution from the lower surface boundary layers 
must be added to complete the solution. Since knowledge of the complementary 
solutions are not needed in the present study, we will not discuss them further. 

In the following section we show how the particular solution, by itself, can be 
employed to correct the external viscid/ inviscid solution for trailing-edge 
interaction effects. A procedure for removing the singularities of second-order 
boundary- layer theory and accounting for normal pressure-gradient effects will also 
be described. 

2.4 COMPOSITE SOLUTION 

In this section we describe a procedure for using the inner trailing-edge 
solution to correct the conventional boundary- layer formulation described in Sec. 
2.2. The composite solution given in Sec. 2-2 is uniformly valid in the inviscid, 
boundary- layer and wake regions but is not valid in the strong-interaction region at 
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trailing edges. A solution, uniformly valid in the trailing edge region, can be 
constructed by adding terms corresponding to the trailing-edge solution minus its 
common part to the "outer" composite solution given in Eqs . (1). This extends the 

domain of validity of the "outer" composite solution to include the trailing-edge 
region. Thus, a composite solution that is uniformly valid in the trailing edge 
region can be written in the following form 

c U=U(s,n; v a ; C V ]] ; r)+ [u(s,n) — u cp (s, H)] +[u(£,tj) — u cp (£ ,tj)] (57a) 

U = V(s, n; v 0 ; [I V E; T) + [v(s, n) — v C p (s, n)] + [v(£, t]) - v cp .(£, rj)] (57b) 

(P = P(s, n; v a ; II V r ; r) + [p(s,B) -p^ (s, n)] + [p(£ , tj) -p cp ( 4 , 7 ])] (57c) 

where u, u , etc. are the contributions from the trailing-edge solution and its 
associated common part. Similar representations are also assumed for the density and 
temperature. As discussed in the previous section, the complete trailing-edge 
solution is the sum of a particular integral and a homogeneous solution. The 
homogeneous solution is a solution of the inviscid irrotational-f low equations that 
is added to the particular integral in order to satisfy boundary conditions in the 
trailing-edge region. The necessity of utilizing homogeneous solutions can be 
avoided if we include only the contributions from the particular integral in the 
expressions for the trailing-edge solution appearing in the last bracket of Eqs. 
(57). With this choice, the matching conditions and the functions v^, [£vj] and V 
which determine the outer inviscid solution, are modified such that the homogeneous 
solution is automatically included in the inviscid outer solution given by the first 
term in Eqs. (57). Thus the trailing-edge terms in Eq. (57) are given by 


v-v cp =e 3/2 v p 
P-P cp = e 3/2 PD 
u-Ucp=e 3/2 S 


(58a) 

(58b) 

(58c) 


where v^ and p^ are determined from the contributions to the particular integral as 
given in Eqs. (53). The term u is obtained from p through the definitions in Eqs. 
(24a), (32b) and (38a). 


The modified matching conditions can be derived from Eqs. (57) by imposing the 
condition that the normal velocity components and v vanish on the airfoil surface 
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(n = n = 0) and the condition that both the full solution ( C U,*U,(P ) and the 
boundary- layer solution (u,v,p) be continuous in the wake. Specializing Eqs . (57a, 
b) to the airfoil surface and wake axis leads to the following expression for the 
composite solution of the normal velocity and pressure on the airfoil surface and 
wake axis: 


U e (s) = V e (s) + € 3/ 2 [v D (| ) + H(-4)J -e 2 v 0 (s) 

(P e (s) = P e (s) + € 3/ 2 j^p p (4 ) + H(§ )1 — e 2 p ff (s) 


(59a), 


(59b) 


where H(£) is the unit step function, and 
° Pe ds 


(60a) 


Pa = ± ^eU e M<5*+^)^(s) (6Qb) 

The plus (minus) sign is for the upper (lower) side of the airfoil and the constants 
A and a^ are given by Eqs. (53c) and (56). 


The condition that the V vanish on the airfoil surface leads to the following 
expression for the corrected source velocity on the airfoil: 

V e (s) = - € V 2 j\ (g ) + ^ Jiff 1 j + e 2 v g (61) 

The requirement that the composite normal velocity and pressure, as determined from 
Eqs. (59), be continuous across the wake axis leads to the following relationships 
for the jumps in source velocity Jv]] and pressure M across the wake axis: 

ttVE = Ve(s) - v;(s) = -e 3 / 2 [Vp(£ + ) - v*(£‘)]+e 2 v w 

(62a) 

ttPB = P e + (s) -P e -(s) = -e 3/2 [p P V ) -Ppin +(^ (1 -a 3 )J-e 2 p w U w c w * 

The n w ,f subscripted quantities are defined by Eqs. (9). 

The terms containing A + and A in the above equations are contributions from the 
common-part terms. Note, there are no common-part contributions to the 
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normal -velocity component terms in the wake or to the pressure distributions on the 
airfoil surface since these common-part contributions vanish to the order of terms 
retained in Eqs. (59), The formula given in Eq. (59) for (P (s) provides composite 
solutions for the pressure on the airfoil surface and wake axis as corrected for 
pressure variations across the boundary layer and wake. 


The matching conditions given in Eqs. (61) and (62) form a sum-type composite 

expansion being given by the sum of a boundary- layer type solution plus the 

trailing-edge solution minus the common part. The added terms from the trailing-edge 

solution cancel the singular contributions appearing in the boundary- layer solution 

(Eqs. 14c, d) resulting in "corrected" viscous matching conditions that are uniformly 

valid in the "outer" boundary- layer and "inner" trailing-edge regions. The- modified 

matching conditions are valid in a formal asymptotic sense for the limit Re ■+ °°. In 

our applications of these equations, the inner trailing edge solution is not carried 

out to the same order as the outer boundary- layer solution. Consequently, in this 

situation, the singular terms in the outer and inner solutions will not be identical, 

3/2 

but will differ by terms that are formally of higher order than 0 (s ). Thus, 
although the inner and outer solutions will formally match term by term if Eqs. (59) 

- (62) are expanded in Re, the singular terms will not exactly match in our method 
where the outer boundary- layer solution is not expanded in a formal asymptotic 
expansion. This would lead to technical problems in numerical solutions based on 
Eqs. (59) - (62). The problem can be avoided and a form suitable for numerical 
computation can be obtained by rearranging part of the matching conditions into a 
multiplicative-type composite solution. Thus, if we consider the boundary- layer and 
trailing edge solutions as the "outer" and "inner" viscous solutions, we can replace 
Eq. (59a) by: 


V e (s) = 


V e (s) — e 2 v 0 (s) 


_ v a.(D J 


V e (s) — e 3/2 av p (£) — e 2 v a (s) 


where 


= V 7 s — s te I K 

v cp 

y - — t — - ) + ) i 

K - J 


(airfoil) 

(wake) 


(63a) 

(63b) 


(64a) 

(64b) 
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3 /2 

The leading order (s ) term in the source velocity given by Eq. (63) is 
continuous across the trailing edge, but the second-order term is discontinuous with 
a small jump across the trailing edge. This jump is a consequence of the 
viscid-inviscid interaction, which acts to smooth the corner formed by the airfoil 
surface and the dividing streamline of the viscous flow leaving the trailing edge 
(see Fig. 8). The value of the constant, a, is determined from the condition that 
the jump in source velocity across the trailing edge is equal to *11(5^) tan &( s te )> 
where 

<U( s tp ) and f5(s te ) are the values of the composite velocity and wake streamline angle 


at the trailing edge. With this condition and Eq. (64), we can rearrange Eq. (63) 
into the form, 


jV e (s) — e v a (s) vl 
V * (S)= \V.(S)-D K(|)- 


V e (s)-€\(S)V I s — s te | K(|) (airfoil) 


e 2 v a (s) (wake) 


where 


K = Ellil+jgzij) 
a t (0) + ia 2 ( 0) 

D = v e(s te ) - e2 v CT (s te ) - < U e (s te )tan/3(s te ) 
where V 


(65a) 

(65b) 


(65c) 

(65d) 


^ (s ) t * ie lilting trailing edge value from the airfoil side as 

determined from Eq. (65a) and v c ( s te ) is the limiting value from the wake side. 


In the above formulae, a multiplicative-type representation is used for the 

velocity on the airfoil surface while a sum-type representation is used in the wake. 

A multiplicative-type composite solution is not possible in the wake because the 

common part vanishes in the wake to the order of the terms retained in Eq. (65). 

However, the above representation is sufficient for our purposes. The square root 

term appearing in Eq. (65a) cancels the corresponding singular term appearing in the 

boundary- layer solution, v fs) for s s, . This leads to a smooth solution near the 

o te 

trailing edge, independent of the order of the terms carried in the inner and outer 

viscous solutions. The parameter, D, is a scaling parameter introduced to assure a 

continuous transition of the inner solution across the trailing edge. The 

representation in Eqs . (65) is formally equivalent to the original sum-type composite 

expansion in Eq. (59a) in the sense that the difference between the two solutions is 
2 

less than 0(e ) . 
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The same considerations lead to the following representation for the airfoil 
surface and wake pressure distribution. Thus, 


<P e (s) - | 


P e (s) -e 3/Z I(£) ttP(s te )D-€ 2 P<J (s) 
P e (s) — II P(s)D 


(airfoil) 

(wake) 


where 


ffPB = e 2 p CT (s) V |s — s te | I(|) 

(TP(s te )D = Pe(s te )-(Pe(s t e) 

1 j~ ffi(£) + jfa(£) j 


I = 


B 1 / 2 6 1 / 2 (l-a 3 ) 


(66a) 

(66b) 

(67a) 

(67b) 

(67c) 

(67d) 


? = M£l±ihl£l 
^(0) + ^ 2 (0) 

The viscous matching conditions can be obtained directly from the above 

representation for the velocity 7J and pressure (P . The condition that the normal 

e 6 

velocity component vanish on the airfoil yields 

V e (0) =e 2 v 0 (s) / I S- Ste | K(4) ( 68 ) 

The requirement that the viscous solution be continuous across the wake axis leads to 

the following expressions for the jump in source velocity [£v]] and pressure JpJ] across 

the wake: ^ __ , _ 

UVD = [D + K"(n-D ITU )]+v w (s) 

(69a) 

IlPB=IlPr-ffPr (69b) 

The matching conditions for the circulation, T, follow from Eq. (69b) and the 
relation 

1IPD + 0(1IPI1 2 ) 


T" = IIU11 = "7V 

ds P w u w 


(70) 


Thus, substitution of Eq. (69b) into Eq. (70) followed by some rearrangement leads to 
the expression 

dr = € 2 c w [Vis — s te I J w ] d/3 (71a) 


where 


and 


t / , r 1 ~i r \ ( ms *) 

Jw L- b1/2 < i — » 3 >J l < 6 *) 1/2 U^+OwA i + ^ + / 

i / fi* + e~ 1 

i‘) l/2 \6* + 6 w /\ 1 + 7?" / J 


( 6 * 


6* = 6* + 6j e Vf =e*+e m 


(71b) 

(71c) 


37 



Equations (68), (69) and (71) constitute the principal results of our viscous theory. 
These equations provide all the boundary conditions needed to determine the 
second-order, outer, inviscid solution. The formulae provide a rational method for 
correcting the matching conditions of conventional interacting boundary layer to 
account for strong-interaction effects near trailing edges of cusped airfoils. 

Before concluding this section we provide alternative formulae for determining 
the pressure on the airfoil and in the wake which are required in the solution of the 
boundary- layer equations. The composite solution for the pressure on the wake axis 
is given by Eq. (66b). Following some straightforward algebraic manipulations, it 
can be expressed in terms of the inviscid quantity P^(s) and a boundary- layer 
function X*, as follows: 

(P e (s) = P e + (s)-Y7^llPB (72a) 


where 


ffPD = P* (s) — P"(s) and \ is defined by 


(6* + e+\ 

"(6 + ) 1/2 (l + ^ + )l 


\6* + e- ) 

L(6-) 1 / 2 (1 + ?-)J 

_01 (!') + ir'/MO J 


The pressure on the airfoil surface can be written in the form 
cP. (s)= p c(s) _ £ 3/2ip (Stt) j [ £;<;; g;*; ] _«»*,.> 

where 

ITP(s te )B = P e (s te ) - c? e (s te ) 


(72b) 


(73a) 


(73b) 
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The last two terms in Eq. (73a) account for the 0(c^ 2 ) and 0(e 2 ) variation in pres- 
sure across the boundary layer. 

2.5 SUMMARY OF THE CORRECTED MATCHING CONDITIONS 

In the theory presented in this report, the inviscid and viscous solutions are 
coupled through viscous matching conditions, just as in the standard interacting 
boundary- layer theory. In the new theory, however, the matching conditions are cor 
rected as described in the previous subsection to account for strong interaction and 
normal pressure-gradient effects in the trailing-edge region. The resulting numer 
ical problem is then very similar to the problem arising in standard interacting 
boundary- layer theory. In both theories the inviscid and boundary- layer equations 
are solved by iteration to obtain self-consistent solutions satisfying the matching 
conditions coupling the two solutions. The new features arising in the present 
method are concerned with corrections to the matching conditions and the need to 
accommodate the wake matching conditions in the numerical solution. The new matching 
conditions are summarized below: 

Source velocity: 


V e ( S )=€ 2 V ff [ 



(airfoil) 

(74a) 

nv u(s) = |d + 

k ( S b#“)- d ' e ( 2 5#“) 

|1 +e2 v w 

(wake) 

(74b) 


Wake circulation: 


r(s) = r 00 +e 2 



c w y[ s — s te J w (s) d jS 


(74c) 
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Pressure distribution: 


(P e (s) = P e (s) — e 3/2 [[p(s te )E I 



± Pe U e K (6* + 0) 


(airfoil) (74d) 


(P e (s)-P e + (s)-^ IIP11 


(wake) (74e) 


where 


M = K " p e" 

The trailing edge corrections used in the above coupling conditions are strictly 
applicable only to airfoils with cusped trailing edges. Although the local strong 
interaction theory can, in principle, be generalized to wedge shaped trailing edges, 
this has not as yet been done. The present theory as formulated in Eqs . (74) leads 
to an equivalent source velocity and flow deflection that has a discontinuity that 
exactly cancels the geometric corner at the trailing edge of an airfoil with a wedge 
shaped trailing edge. This is exactly the same behavior as in standard interacting 
boundary layer theory. This discontinuity in source velocity across the trailing 
edge requires us to modify the determination of the constant D defined by Eq. (65d) . 
For airfoils with included trailing edge angles, 0^, the equation for the scaling 
parameter, D, should be changed to, 

D* = V*(s to ) - e 2 v£(s ta) - ‘Ws te ) tan [0 (s te ) ±0 to / 2 ] (74f ) 

where the last term accounts for the jump in source velocity due to the nonzero 
trailing edge angle. 
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3. OUTER (POTENTIAL FLOW) SOLUTION 


In this section we outline the procedures used to solve the outer inviscid 
problem, assuming the boundary- layer quantities appearing in the matching conditions 
are known from a previous iteration. The solution in the outer inviscid region is 
based on a potential-flow approximation which should be adequate for weak shock 
waves. We employ Jameson’s (Refs. 6 & 7) hybrid, relaxation/direct-solver scheme for 
solving the full-potential equation in conservative form. The particular version we 
employ is described in Ref. 7. For completeness, we briefly review Jameson’s scheme 
and indicate where modifications are necessary to accommodate the viscous matching 
conditions . 

3.1 GOVERNING EQUATIONS AND BOUNDARY CONDITIONS 

Jameson’s method is carried out in a computational plane that is developed from 
the conformal transformation of the region exterior to the airfoil to the interior of 
the unit circle. A polar coordinate system, (r, w) is used in the circle plane to 
generate a desirable grid system for the finite-difference approximation. The 
potential-flow problem to be solved in the circle plane is indicated in Fig. 7. The 
conservation equation for the velocity potential in the circle plane can be written 
in the form 


^(KQ„) + r i i(EQ f ) = 0 


where 


Qw 


9 $ 

0U’ ’ 



(75) 


(76) 


and where, in the notation of the previous section, R is the density in the outer 
(inviscid region) and Q_ and are the mapped velocity components in the r and u 
directions, respectively. If3C/r is the modulus of the transformation of the airfoil 
into the exterior of the circle, the physical velocity components in the F and u 
directions may be written as 
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AIRFOIL ("= 1) 



Figure 7 Computational Domain for the Inviscid Problem 






(77) 


The density, R, is determined from speed of sound, a, and the energy equation in 
isentropic flow, as follows 


a 2 = 




*0 


(Qf + q^) 


(78a) 


R = (Mia 2 ) 1/r_1 


(78b) 


where a^ is the stagnation speed of sound. The uniform free stream and the net mass 
flow from the boundary layer introduce singularities at the origin, r = 0 , in the 
circle plane. These are removed by introducing a reduced potential, G, defined by 


~ _ cos(a' + a) 

G = <p — r 

r 


-e 2 cj d ln r + E(co + a) 


(79a) 


where 2 ttE is the circulation at infinity and is the source flow at infinity 
introduced by the boundary layer and wake. Thus, 


E = lim 

r** 0 


r(Cr) 

2tt 


(79b) 
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(79c) 


e 2 (r d = lim 

r -*0 


27tB«, 


where 


Boo = V 1 — M 4 


(79d) 


The transformation leads to a reduced potential that is single-valued and bounded at 
infinity. The corresponding far-field boundary condition is 


lim G= E{ (to + a) — tan’ 1 [ B*, tan(o: + a)] } — — ^ In [ 1 — sin 2 (o: + a)] 
r -0 1 


(80) 


The last term in Eq. (80) arises from the viscous -induced source flow at infinity. 
The mapped velocity components and can be written in terms of the reduced 
potential as 

3G sin(a? + Qf) 

Q ““9co ~ r (81a) 

_ aG_cos(^+a) +e2CT ( 81b) 

^ r 3r r d 

The boundary condition on the airfoil surface (r = 1) is 

= cos(a' + a) — € 2 cr d + 3CV e (82) 

where is the surface-source velocity defined in Eqs . (68). 

To be consistent with the theoretical analysis, the wake boundary conditions 
should be imposed along a reference curve that is close to the wake, say within 0 (s) 
from the wake centerline. We stress the fact that the reference wake axis need not 
be a streamline, but that it must only meet the requirement that it be sufficiently 
close to the wake centerline. In the present study we align the reference wake axis 
with the radial cut w = 0 in the circle plane as sketched in Fig. 8. Both the "cut" 
and the inviscid streamline are tangent to the airfoil bisector at the trailing edge. 
With this choice, the wake axis is tangent to the trailing streamline of the inviscid 
solution and will, therefore, be positioned close to the wake near the trailing edge. 
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Figure 8 Wake Reference Axis 

Since the wake matching terms decrease rapidly away from the trailing edge and 
are quite small at distances greater than the present chord, the approximation should 
be more than adequate. Thus, the matching conditions in the wake along w = 0 can be 
written in the form 



E+ sin(u> + a) + 3C[[V]] 


(83a) 


G + — G" =Ar(r) = Ar(r)-r(0) 


(83b) 


where [V] is the jump in source velocity given by Eq. (69a) and T(r) is the 
circulation distribution along the wake which is determined from the integration of 
Eq. (17a) subject to the boundary condition at r = 0 

r(0) = 2irE ( 84 ) 

The circulation at infinity (r = 0) and the constant E are determined from the Kutta 
condition requiring the velocity at the trailing edge to be bounded. This is 
satisfied by requiring Q y to be zero at the trailing edge. Therefore, from Eq. (81) 
we arrive at the following condition for E 

9G 

E = sinot + — atr =1, 0 (85) 

The determination of Y (r) from Eq. (71a) requires knowledge of the wake angle, 

In the present work B is defined to be equal to the angle between the extension of 
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the airfoil chord and the composite streamline leaving the trailing edge (streamline 
1 in Fig. 8). If this streamline is assumed to be close to the position of the cut, 
then a sufficiently accurate approximation is 

0 = + O (ycut'yo) (86) 

where 3 r is the angle between the radial line w = 0 and the extension of the airfoil 
chord and 3 q is the local flow angle at the "cut" as depicted in Fig. 8. The local 
flow angle, 3 q 5 can be determined from the composite solutions for the two velocity 
components V e ,TL e at the cut, w = 0; thus, by definition 

= tan _ 1 CU e >/%) (87) 

where is evaluated from Eq. (65b) using the modified expression for the scaling 

e 

parameter D given in Eq. (74f) 

^ e (S) = V e (s)-DK({)-e\(s) (88) 

The error term, (y ut ~ y Q ) appearing in Eq. (86) is the distance between the cut and 
the wake centerline. For inviscid flow, the terms c^, and [[v]] are zero, the 
circulation, T, is equal to 2ttE and the above formulation reduces exactly to 
Jameson’s original fully-conservative scheme for the full potential equation. 

3.2 NUMERICAL SOLUTION 

The outer inviscid problem formulated above is solved by Jameson’s (Refs. 6 
and 7) fully-conservative relaxation method. The method employs a rotated-dif ference 
scheme and a convergence acceleration technique based on a combined SLOR-Poisson 
direct solver. Convergence is improved by carrying out the computations on a 
sequence of three meshes. Options for a standard non-conservative formulation are 
also programmed into the method. The basic numerical method is fully described in 
the original references (Refs. 6 and 7). Here, we describe only the changes needed 
to accommodate the viscous terms in the boundary conditions. 

The computational grid employed in the calculations is indicated in Fig. 9 
together with a sketch illustrating the grid distribution in the physical plane. The 
flow field in the computational plane is contained between the airfoil surface 7=1, 
the "point" at infinity f = 0, and the two sides of the "cut" w = 0 , w = 2 tt . The 
computational mesh consists of interior points, image points used to satisfy the 
airfoil boundary condition and overlap points employed to satisfy the wake matching 
conditions. The indices i, j are used to label grid points on the w, r axes, 
respectively. The present method overlaps the grid on both the upper and lower sides 
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of the cut, in contrast to Jameson* s original method which uses overlap points only 
on the lower side of the cut. 


The body boundary condition is imposed by inserting appropriate values of mass 
flux normal to the airfoil surface, RQ_, into the image points. The values at the 
image points (i, J+l) are obtained by setting the mass flux at the surface, j = J, 
equal to the average of the points (J + 1) and (J - 1). Thus, 

(RQr) lf J + l - - (RQ?) lf j-l -2 OCR^e (89) 


The wake conditions are imposed by setting appropriate values of the potential 

into the overlap, i = 2, and !t cut M , rows i = N. The solution in the overlap row, 

i = N + 1, is used in the derivation of the wake conditions but is not actually 

employed in the numerical solution. The arrangement of the grid points across the 

cut is indicated in Figs. 9 and 10. The reflection rules for the determination of 

the potentials ^ and ^ can be found from the finite-difference expressions for 

the jumps in G, G and G across the cut. Thus, from the definition of T we have 
w ww 5 



OVERLAPPED 

RAYS 
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G N,j “ G 3,j 


where AI\ is the value of Ar at r = r.. By using centered-dif ference expressions for 
on both sides of the cut in Eq. (83a), we can obtain the following expression for 
the mass -flux discontinuity across the cut 

(G n+ I,j “ g n-i, j) - ( G 4,j ~ G 2,j) = “ 2Acr3C[[Vl] (91) 

Similarly, an expression for the jump in G ■ across the cut can be obtained in the 


form 

(G n 

where 

(Tg 


IL tow 

F=l 

3 

3 

L=l 

is 

terms 

of G. 


G 2,j = Gn-i,j “ Ar j + Aoj 3 C|[v ]1 - 5 Aw 2 llG UC0 I] + 0( Aw®) 
G n+1>j = G N . lfJ +AT i AcaCffVl + iAco 2 I[G ww B + 0(Aa< 3 ) 


(93a) 

(93b) 


The solution is carried out by sweeping the field along columns i = constant 

from the leading edge towards the trailing edge, first on the lower surface and then 

on the upper surface. In carrying out the sweep on the lower surface, the value of 

G^ j is used from Eq. (93a) with G^_^ j and r evaluated from the previous sweep. The 

sweep on the upper surface employs the potential G M . as determined from Eq. (90) as 

M 5 J 

a boundary value in the solution. The circulation is obtained from an integration of 
Eq. (71a) using a simple central difference approximation to yield 


Fj+l “ Fj + ” J| c w l s “ s te I ^w)r=rj +1 + (c w VTs ®teT^w)r=rJ ^J+l 


(94a) 

with 

r 2 =27rE (94b) 

The value of E is determined from the Kutta condition, expressed in finite-difference 
form as 


Ga t “ Go 


— sina; + 0(Aor) 


2 Ao? ' — 7 (95) 

The second-derivative term, [°«j , could be determined from the differential 
equation evaluated on each side of the cut. However, we follow a simpler procedure 
and evaluate this term numerically from central difference formulae centered one grid 
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point on either side of the cut. This provides a first-order accurate estimate for 
I G uJ which is sufficient to evaluate the second-order terms in Eq. (93). Thus, we 
can write 

= (G wcj )n-i,j “ (G W a))4 f j +0(Aw) (96a) 

and hence 

IIG WW 11 = [ (Gnj - 2GN. ltJ +G N _ 2fJ ) - (G 5|J - 2G 4fj + G 3>j )] / Aw 2 + 0(Aw) (96b) 


The flow angle appearing in the equation for the circulation, Eq. (94a), is 
determined from the velocities on the cut through Eqs . (86) - (88). To minimize 

numerical error we compute *1)^ as the average of *U~(s) ° n u PP er anc * l° wer 

surfaces of the wake. Thus U e (s) in Eq. (87) is computed from the following 
expression 

V e (s)=h[vt(s)+v- e ( s )] 

+ (s) and V"e a PP ear i n S in Eqs. (88) and (97) are the velocity components 

where V e ^ ^ 

normal to the cut (in Jameson's notation; the q^ component defined by Eq. (77)). 
These can be expressed in terras of the reduced potential by a central-difference 
expression of G centered about the column i = N and i = 3, respectively. Thus, 


v * +(s)=_ ^ 
Ve ‘ (s) =~x 



(98a) 

(98b) 


where the potentials G^ + ^ and G^ ^ appearing in the above equations are evaluated 
from Eq. (93). The evaluation of the composite expression for radial velocity 
component ,qi , required in Eq. (87) is described in Section 4. 


The formula for determining the velocities from the reduced potential, Eq. (77) 
are indeterminate at the trailing edge since the metric, 3C, vanishes at this point. 
Hence, special formulae are needed to compute both the velocity and potential at thi: 
point. In the original inviscid method (Refs. 6 and 7) the velocity at the trailing 
edge was determined from a simple extrapolation of the velocity from neighboring 
points on the airfoil side of the trailing edge. In the present work we employ a 
simple extrapolation from the wake side. Since the velocity variations are more 
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gradual in the wake, this change provides an improved determination of the velocity 
at the trailing edge. A special formula is also needed to evaluate the potential at 
the trailing edge since the transformed potential-flow equation is also indeterminate 
at this point. We employ the simple interpolation formula used in the inviscid 
solution in Refs. 6 and 7 based on the grid points indicated in the Fig. 11. Thus, 
with this scheme the trailing-edge potential, G., . is determined from the formula 

^ 5 J 

G 4,J “ 2G 3 f J + G 2,J G 3 f J-l~ 2G 3,J +G 3 f J+l _ „ 

+ “ u (99a) 

A somewhat more complicated procedure based on the use of an approximate 

Prandtl-Glauert equation in the physical plane near the trailing edge is employed in 

Ref. (17). 

The difference equations are solved for fixed estimates of the boundary- layer 
parameters V^ y Y , and Jj^vj] . The equations are solved by the fast iterative scheme 
described in Refs. 6 and 7 which is based on alternating a sequence of SLOR and 
Poisson iterative methods. In each work cycle, we first perform a Poisson step, 
followed by a specified number of relaxation steps. During the computation we 
monitor the ratio of the maximum residuals at the start and end of a work cycle. 
When this "residual ratio" is reduced to a prescribed level, the boundary- layer 
calculations are repeated and the boundary- layer parameters appearing in the matching 
conditions are updated. The inviscid computations are then repeated. This 
alternating sequence of inviscid and boundary layer computations is continued until a 
set of convergence criteria are satisfied. The criteria employed are based on the 



“CUT" 
(cj = 0) 


Figure 1 1 Interpolation Points for the Trailing Edge Potential 
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maximum residual in the inviscid computation and the change in lift coefficient (a 
prescribed) or angle of attack (C T prescribed) between successive work cycles. We 
also provide for a termination of the computation based on a maximum number of 
iterations. The calculations are carried out on a sequence of three meshes. The 
solution on the crude grid is initiated using built-in estimates for the 
boundary- layer parameters and reduced potential. The solutions on succeeding grids 
are initiated from interpolated boundary- layer parameters and reduced potentials. 
The solution on each grid starts with a solution of the inviscid equations. 
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4. BOUNDARY LAYER SOLUTION 


The boundary- layer parameters appearing in the matching conditions are evaluated 
from solution to the boundary- layer equations with a known pressure distribution 
determined from the previous iterate. Since the present method accounts for pressure 
variations across the boundary layer and wake, there is a certain degree of arbi- 
trariness in the choice of the pressure distribution to be used in the boundary layer 
computations. We could, for example, use the pressure distribution at either the 
edge of the boundary layer or at the surface or, could use an appropriate average 

across the boundary layer. Since the pressure variations across the boundary layer 

3/2 

and wake are small, being at most of 0 (s ) near the trailing edge, the various 
choices lead to equivalent solutions, to lowest order in £. In the present work we 
use the composite solution for the pressure, (p^(s), on the airfoil surface, Eq. (73), 
and cut, Eq. (72a), neglecting the small second-order term P^(s) appearing in Eq. 
(73.) In the solution of the boundary layer equations it is assumed the pressure is 
constant across the boundary layer and is given by(p^(s). 

The edge velocity, U^ ^(s), and density, ^(s), appearing in the boundary 
layer equations and matching conditions can be evaluated from the exact isentropic 
relations . 


U e 2 rW (s) = l + 


[yMj(P.(8)l y ^-l 
(y- 1)M i/2 


(100a) 


Pe,w(s) = [yMi<P e (s)] 1/r 


(100b) 


In the evaluation of the flow angle (5 in Eq. (87), we approximate 'll (s) by 

O V 

U^(s) as determined from Eq. (100a) which is a consistent approximation to lowest 
order in z . 


In the present work both the laminar and turbulent boundary- layer equations are 
solved using simple integral methods. In this section we briefly describe the 
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methods used to solve the boundary- layer equations and the procedures employed for 
transition and leading edge separation. 

4.1 LAMINAR BOUNDARY LAYER, LEADING-EDGE SEPARATION AND TRANSITION 

The laminar boundary- layer equations are solved for 6*, 8, h and by a 

compressible version (Ref. 9) of Thwaites 1 integral method (Ref. 61) for incompress- 
ible flow. The laminar equations are integrated from the forward stagnation point to 
a transition point whose location is either assigned or determined by a specified 
transition criteria. Three semi-empirical, transition prediction methods have been 
programmed into our method as user selected options. These include transition crite- 
ria based on Crabtree ! s correlation (Ref. 62), Michel ! s correlation (Ref. 63), and 
the correlation in Stevens, Goradia, and Braden (Ref. 64). We also test the laminar 
solution for leading-edge separation according to the criterion in Ref. 64. If 
leading-edge separation is present, a diagnostic is written indicating that sepa- 
ration has occurred and whether it is of the long or short bubble type. If transi- 
tion is not predicted before laminar separation, transition to turbulent flow is 
assumed to occur at the laminar separation point. The user may choose both an 
assigned transition location and one of the above transition criteria. In this case 
the transition point used in the calculations is the most upstream of 1) the assigned 
location, 2) the position predicted by the transition criteria or 3) the laminar 
separation point. It is, of course, well known that these and other existing methods 
of predicting transition are not reliable. They are included in the present method 
merely as a guide to be used in the absence of any information regarding the position 
of transition. 

The turbulent boundary- layer calculation is initiated at the transition point, 
generally, under the assumption that the momentum thickness is continuous and the 
shape factor is discontinuous across transition. The initial value of the shape 
factor in the turbulent flow is determined under the assumption that the increment in 
the "incompressible” shape factor, Ti, is given by 

Ah t =1.1 < 101 ) 

where, in turbulent flow 


h = (h.+ l)(l+ r -^f) -1 


( 102 ) 
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where r is the recovery factor and is the local Mach number at the boundary layer 
edge. A more elaborate correlation in Schlichting (Ref. 65) based on the assumption 
that Ah t is a universal function of the momentum thickness Reynolds number at transi- 
tion, can be used in place of Eq. (101). However, the turbulent boundary layer down- 
stream of transition is not sensitive to the initial conditions at transition and the 
simple relationship given in Eq. (101) is adequate. In our method, we also allow for 
a jump in momentum thickness, AG^., across the transition point. This jump can be 
used to simulate an increase in momentum thickness due to a roughness strip. Thus, 
the initial values of the reduced shape factor and momentum thickness employed in the 
calculation of the turbulent boundary layer are written in the form 


h - hi am i nar (co t ) + Ah t (103a) 

0 = 0 lamla Ju t ) + Ae t (103b) 


where u = is the circle plane angle of the transition point. 


The integrals appearing in the laminar solution are evaluated by a simple 
trapezoidal rule and the surface source velocity at points where the flow is laminar 
is evaluated from the following central-difference approximation: 


V e (co t ) = 


( 1 ) 

f (p e u e 6 * ) i .1 - (p e u e 6 * ) ! ., 1 

VPeK ) 

iL 


4.2 TURBULENT BOUNDARY LAYER SOLUTION 


(104) 


The solution of the turbulent boundary layer and wake are carried out by the 
lag-entrainment method of Green et al. (Ref. 10). This is an entrainment -type inte- 
gral method that includes an approximate treatment of the turbulent energy equation. 
The method has been demonstrated to provide accurate solutions for airfoil-type flows 
with large adverse pressure gradients. Such flows arise, for example, near shock 
waves and trailing edges of rear- loaded supercritical airfoils. 


The basic method employs the momentum- integral equation, a shape-factor equa- 
tion, and a differential equation for the entrainment function. To avoid interpolat- 
ing between the inviscid and boundary- layer solutions, the latter are integrated in 
the circle plane using the same grid points employed in the inviscid solution. The 
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basic equations for the momentum thickness, reduced shape factor, and entrainment 
function and are written in the form 


i_ de 
3C d£ 


= ic t — (h + 2 — M|) 


e 

U e 3C 


diL 

dl 


1 dh 1 dh 
3C d£ " e dh t 



e du 0 

(h+1) IL d£ 


(105) 


(106) 


1 dc E _ 1 F 

0C d£ e e 


(107) 


where Z is equal to the polar angle w on the airfoil surface and to the radius, r, in 
the wake (We recall that 3C is the modulus of the transformation.)* The shape fac- 
tor, h, is determined from h by Eq. (102) and the displacement thickness is given by 


6* = h0 


(108) 


The velocity-profile function, h^, is given in terms of h by 

h, = 3.15 + ~ .01 (h-1.) 2 

(h-1) 


(109) 


The skin friction is computed from the following empirical relations 


, -e r _ ,_1 1 
" f f0 [_10h-4 2 J 


( 110 ) 


where 


h 0 = [ 1 + 6. 55 V| Cf (1+.04M^))" 1/2 

1 f 0.01013 i 

Cf0 F c [lo glo (F R R 5 )-1.02 • 0075 J 

F c = VI + . 2M| } F R = 1 . 0 + 0. 056 

and Rg is the local Reynolds number based on momentum thickness. The function 
appearing in the entrainment equation is an empirical function defined in Ref. 10 

that depends on boundary- layer parameters, the external-velocity gradient dU g 

Ue ds * 

the local external Mach number M . This function also involves parameters that per- 

e 

mit a rough treatment of the effects of longitudinal curvature and mean dilatation 


(111a) 

(111b) 

(111c) 
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on the turbulence. The same model is used to solve for both the flow in the boundary 
layer and the wake. The wake is analyzed as if it were two separate symmetric half 
wakes . 


The integration of Eqs . (105) - (107) is initiated at the transition point using 
values for h and 0 determined from Eqs. (103). The initial value of the entrainment 
function, C^, is determined from an equilibrium value defined in Ref. 10. For fur- 
ther details of the model see Ref. 10. 


The differential equations are integrated using a standard variable-step Runge- 
Kutta method. The integrations are carried out on a basic grid formed by the w nodal 
values used in the inviscid solution. The integration scheme allows repeated grid 
halving, subject to an error measure, to maintain accuracy in regions of rapid varia- 
tion as occur near shock waves and the trailing edge. In this way the solution for 
the boundary- layer parameters are determined directly at the nodal points of the 
inviscid flow without interpolation. 


A useful expression for the source velocity, v^, in turbulent flow can be ob- 
tained from the momentum- integral and shape-factor equations. Thus, employing Eqs. 
(105) and (106) in the definition of v^ we arrive at the following expression for v^ 
as a function of the gradient of the edge velocity U^: 

dU e 

V(T = Cl_C2 “d7 (112) 


where 


c i = 


Y hc f + (c E -|h 1 c f )(c r +l) 



c 2 


= 0(h+l) 



(e E — 1) (Cp + 1) 
(c r + 1) 


- 2 < C " +1) 


h dh 1 

h 'dhTj 


and 


c r = 1 + r(y — 1)M| , c 0 =1 + (Y-1)M 2 

dh (h-1) 2 

dh t 1.72 + 0.02(H-1) J 


(113a) 

(113b) 


(113c) 

(113d) 
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In the numerical method, the velocity derivative appearing in Eq. (112) is approx- 
imated by a simple centered difference in the circle plane. It is of some interest 
to note that in the limit Re <», h 1 and approaches the following limit (for 
recovery factor r = 1) 


r 1 / 

f i\„i 

dUel 

Ls c ° Cf 

C 0 + ) TT 

K 0 Co/ U e 

ds J 


(114) 


The boundary layer thickness, 6, and Coles 1 wake parameter, tt , are determined by 
matching the skin friction and displacement thickness corresponding to Coles 1 com- 
pressible law-of -the-wake (Refs. 66 & 67) to the solution of the lag equations just 
upstream of the trailing edge. 


The lag-entrainment method was not designed to deal with separation of the 
boundary layer. In many cases of interest, small separation zones will appear in 
airfoils flows, without significant influence on the surface-pressure distribution or 
other section characteristics. For example, small separation regions can occur at 
the foot of the shock wave or near the trailing edge without exerting a strong influ- 
ence on the outer inviscid flow. In the case of shock- induced separation, the shape 
factor is observed to undergo a very rapid increase as the shock wave is approached, 
followed by an even steeper decrease behind the shock wave. In these cases, the lag- 
entrainment method will predict an increase in the shape factor to unrealistically 
large values (e.g., h 15), completely outside the range of data correlations on 
which the method is founded. The method will also predict unrealistic values of the 
skin friction in the separated zone. In order to enable the present method to func- 
tion in these cases, we have set arbitrary bounds on both the skin friction and shape 
factor when the flow is separated. Flow separation is predicted in the theory when 
the skin friction computed by the lag -entrainment method is zero. If separation is 
predicted, we set an upper bound on h given by h . Then, if h computed according 

to Eq. (106) is greater than h , it is set equal to h . The value of the upper 

max ^ max 

bound can be specified as part of the input; typically, we take h =4.0. For the 

max 

skin friction, we follow Hunter and Reeves (Ref. 68) in their !1 wake like 11 model of 
separation and set c^ = 0 if it is predicted to be less than zero in the lag-entrain- 
ment solution. This approximation of the skin friction in separation zones can be 
expected to have only a small influence on the solution since c^ is generally small 
in the slender separation bubbles that occur on airfoils. The overall accuracy of 
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this approach to separation is not expected to be very high but it should provide 
useful estimates of incipient separation at the shock wave or trailing edge. 


58 



5 . RESULTS 


The theoretical method described in previous sections has been developed into a 
computer code GRUMFOIL to solve the viscous transonic flow over airfoils. A complete 
description of the most recent version of the code is given in Ref. 69 along with 
instructions for its operation. In the version of the code used for the present 
computations, we allow for the choice of the constant, which controls the 

accuracy of the difference approximation in the supersonic zone and for the choice of 
the parameter, Q , which determines the Mach number of the switch from a central to a 
backward difference. The difference formulae are second (first) -order accurate at 
supersonic points for c^ = 1 (0) independent of the value of Q^. In this section 
comparisons of theoretical solutions for various values of c^ and Q are presented to 
demonstrate the effect of the difference approximation on the solution. We also 
present solutions that illustrate the effect of switching from a F-C or N-C 
difference scheme. In the program we also provide an option for deleting the 
trailing-edge corrections and for selectively dropping each of the terms appearing in 
the viscous matching conditions. These options are used to generate results which 
illustrate the effect of wake curvature, wake thickness, and trailing-edge correction 
on the solution. 

Although the basic theory used in the present work formally applies only to 
closed airfoils with cusped trailing edges, the code can be applied to more general 
airfoils not satisfying these conditions. Open trailing-edge airfoils are modeled by 
continuing the airfoil surface downstream with a semi-infinite streamtube. The 
region outside the airfoil plus extensions is then mapped to the region outside a 
unit circle with the two edges of the wake streamtube mapped to a single curve in the 
circle plane. This procedure leads to a constant -thickness airfoil extension and to 
a solution that is continuous across the wake streamtube. For an open airfoil, we 
must also prescribe the base pressure and determine the base-drag contribution to the 
total drag. In the present work we assume the base pressure is equal to the pressure 
at the trailing edge as given by the viscous solution. This approximation is valid 
if the base thickness is less than the boundary- layer momentum thickness at the \ 
trailing edge. The base drag contribution is usually negative, leading to a 
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reduction in drag. The contribution is sensitive to the accuracy of the predicted 
trailing edge pressure and can be significant. Unfortunately, the base drag 
contribution has been missed in recent investigations (Refs. 12-15). 

In this section we also present comparisons of theoretical solutions with 

experimental data for a series of airfoils tested in various wind tunnels. The 

influence of wall-interference effects on the wind-tunnel data introduces 

considerable uncertainty into the interpretation of these comparisons. For small 

enough models, the main effect of wall interference can be taken into account through 

standard downwash and blockage corrections to the experimental values of the 

incidence and free-stream Mach numbers. Unfortunately, reliable estimates for these 

corrections are not available for most transonic wind tunnels, including all those 

considered in the present study. Therefore, in the present investigation we avoid 

interpretational difficulties associated with the angle of attack by comparing theory 

and experiment at the same (measured) lift coefficient. In these calculations the 

lift coefficient is prescribed and the angle of attack is adjusted during the 

relaxation process to satisfy the Kutta-condition at the trailing edge. If reliable 

incidence corrections were available for a set of data, the difference between the 

experimental and theoretical incidence could be used to judge the adequacy of the 

theory. However, because of the absence of reliable incidence corrections, we have 

not attempted such comparisons in the present study. The problem with blockage 

corrections also remains. The main effect of wall-induced Mach number corrections is 

to alter the position of shock waves on the airfoil surface. Mach number corrections 

as small as M = 0.001 can produce noticeable shifts in shock-wave location. When 

available, we followed the recommendations of tunnel operators in applying Mach 

number corrections to the data. Unfortunately, accurate blockage corrections to the 

Mach number were not generally available for most of the wind-tunnel data considered 

in the present study. Therefore, in carrying out the comparisons presented later in 

this section, we employed Mach number shifts for each tunnel that produced the best 

overall agreement between the theoretical and experimental shock locations . Because 

of the uncertainty in the free-stream Mach number, the conclusions drawn from the 

present study must be regarded as tentative. More definitive evaluations of the 

theory must await the availability of interference-free data or accurate 

wall-interference corrections. All theoretical calculations presented in this report 

were carried out on a sequence of three (N x N ) grids consisting of (40 x 8) , (80 x 

r w 
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16) and (160 x 32) points. Except where noted, all calculations were carried out 

with the fully-conservative, first-order accurate ( c 2 - 0, Q =1) version of the 
method. 

5 • 1 effect OF INVIS CID FINITE-DIFFERENCE PARAMETERS ON SOLUTION ACCURACY 

The formal accuracy of the windward finite-difference approximation employed in 

and near the supersonic zone is controlled by the parameter c 2> with the solution 

being first-order accurate for c 2 = 0 and second-order for c 2 = 1. Here we present 

typical results which illustrate the effect of c 2 on the solution. In general, 

convergence slows and the solution becomes more difficult to obtain as c 2 is 

increased from zero to one. Fortunately, convergence is facilitated by reducing the 

value of the parameter, Q c> which controls the points at which the numerical scheme 

is switched from central to backward differences. In the basic setup, Q = 1 and the 

switch is made at the sonic line (M = Q c = 1). Reducing Q c below one tends to smooth 

the solution near shock waves and to speedup convergence. This smoothing effect of 

Q c is especially helpful for values of c 2 close to one. We have been able to 

routinely obtain fully second-order results (c = 1) with Q =0.9. We should stress 

z c 

t at the f inite-dif f erence approximation is formally second-order accurate for c = 1 

independent of the value of 0 . 

c 

To illustrate the effect of c 2 on inviscid solutions we carried out 
fully-conservative computations on an RAE 2822 airfoil at M = 0.725 for three values 
of C L with various combinations of c £ and Q c . This resulted in supercritical flow 

with shock waves in all cases. The results for Cp, a, and shock position, X are 

summarized in Table 1. S 

We first note, the choice of Q £ has only a minimal effect on the solutions for 
fixed values of Second, and more importantly, we note that the solutions for all 

three values of C L> are only weakly dependent on the value of c 2 - The drag 
coefficient is seen to increase by only two counts in going from a first-order (c 2 = 
0) to a second-order (c 2 = 1) scheme. The shock position is uneffected and the angle 
of attack is only slightly increased (Aa .01) as c 2 is increased from zero to one. 

The results also indicate that the solution for c £ = 0.8 and 0.9 are virtually 

identical to the second-order (c 2 = 1) solution. The pressure distribution on the 
airfoil surface for the C L = 0.95 cases given in Fig. 12. The two solutions are 
indistinguishable in this plot. These results indicate that the first-order accurate 
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TABLE lo EFFECT OF c 2 AND Q c ON SOLUTIONS FOR INVISCID FLOW 
OVER A RATE 2822 AIRFOIL AT M = 0.725 

CO 


C L 

C 2 

Q c 

C D 

a 

X 

s 

0.9 

0 

1 

0058 

1.844 

.65 

0.9 

0 

.90 

0058 

1.853 

.65 

0.9 

.8 

.90 

0060 

1.830 

.65 

0.9 

08 

.95 

0060 

1.829 

.65 

0.9 

.9 

.90 

0060 

1.827 

.65 

0.9 

.9 

.95 

0060 

1.827 

. 65 

0.9 

1.0 

.90 

0060 

1.824 

.65 

.925 

0 

1 

0069 

1.904 

.66 

.925 

.90 

.90 

0071 

1.885 

.66 

.925 

1.0 

.90 

0071 

1 o 881 

i 

.66 

.950 

0 

1 

0082 

1 

1.960 

.66 

.950 

.9 

.9 

0084 

1.938 

.66 


version of the F-C scheme employed in the present method is surprisingly accurate 
suggesting that the first-order truncation error terms must be relatively small. Our 
experience with the present method is at variance with the results of Collyer and 
Lock (Ref. 18), obtained for a similar case with their partially-conservative P-C 
scheme. Their results showed a much larger effect of c^, with the drag increasing by 
about twenty counts as c^ increased from zero to one. The major discrepancy is 
between the two first-order solutions, with Collyer and Lock's solutions producing 
much lower drags. The reason for the poor performance of the Collyer and Lock 
first-order P-C method has not been ascertained. 

The above results indicate that the first-order accurate version of the present 
method should be adequate for most purposes. However, there are some cases involving 
weak shock waves for which variations of c^ produce more noticeable effects. In 
these cases weak shock waves near the leading edge were not adequately resolved with 
the first-order accurate version of the method. We found improved resolution could 
be obtained by either refining the mesh or by increasing c^ toward one. A typical 
result illustrating this effect in a full viscous solution is given in Fig. 13. In 
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this figure first and second-order accurate solutions for the surface-pressure 
distributions are compared for M = 0.725 and C T = 0.521. This case involves two 

oo jj 

shock waves on the upper surface of the airfoil. The results indicate that the 
forward, weak shock wave is not adequately resolved by the first-order scheme on the 
(160 x 32) grid employed in the computation, but that good resolution is achieved 
with the second-order scheme. Note that, as is generally the case, the drag and 
angle of attack are only slightly affected by the choice of c^. 

The dotted line extending from the airfoil trailing edge in Figure 13 (and in 
all others plots to follow) is the transformed "cut" (w = 0) location in the 
physical plane. The nearby solid curve is an approximation to the composite solution 
for the streamline from the trailing edge as determined from an itegration of the 
flow angle, f$ “ + Pq, along the "cut" . The closeness of these two curves near the 

trailing edge is supportive of the approximations used in the evaluation of the 
wake-curvature terms. 

The first-order scheme is less expensive to use and has been found to be 
generally adequate for most cases. Therefore, it is recommended for general use and 
has been employed in all solutions presented in this report. The second-order scheme 
is reserved for certain special cases requiring greater accuracy to resolve fine 
details of the flow field as in Fig. 13. 

5.2 EFFECTS OF THE WAKE AND TRAILING EDGE INTERACTION 

The computer code has been organized so that individual terms appearing in the 
viscous matching conditions can be selectively dropped from the computations. This 
option has been employed to generate a series of solutions that illustrate the 
importance of the individual terms in the matching conditions. Three cases have been 
carried out, corresponding to subcritical and supercritical flows over an RAE 2822 
airfoil and to supercritical flow over a more heavily rear- loaded supercritical 
airfoil developed by the NASA (LaRC). For each case, a series of solutions have been 
obtained with each of the terms appearing in the matching conditions selectively 
dropped. Solutions for the drag, lift, and trailing edge pressure are summarized in 
Tables 2A and 2B. In these tables is the total drag, determined from integration 

of the pressure and skin friction over the airfoil surface, and C is the profile 

Jjoo 

drag determined from the wake momentum thickness far downstream (i.e., CL. = 20 ) . 

Jjoo oo 

P r °fil e drag is equal to the total drag less the wave drag and is due solely to 
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-i .600 r 


-1 .200 b 


C 2 =1 C D = 0.0076 a =0.140 
c 2 = 0 C D = 0.0075 a = 0.145 


-0.800 


C p -0.400 


0.000 


0.400 


0.800 


1.200 



momentum losses in the boundary layers For subcritical flow, is equal to the 

total drag of airfoil and should be equal to CL,-, as determined by surface 

UD 

integration. In general, however, these two evaluations of the drag will not be 
equal with the differences due both to, numerical errors in the solution of the 
governing equations and to approximations in the formulation of the viscous effects. 
The momentum drag (i.e., 28 ) seems to be relatively insensitive to the details of 

oo 

the numerical method and the formulation of the viscous theory, and is thought to be 
the more accurate prediction of drag in subcritical flows. In these cases the 
difference between profile and integrated drag is a useful measure of the overall 
accuracy of the theoretical model. 
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The results in the tables include solutions for (A) the inviscid theory, (B) the 
full viscous theory, (C) the full viscous theory less the trailing-edge corrections, 
(D) less the wake-curvature term, and (E) less the wake-thickness terms. In the 
latter solution the wake is modeled as a constant thickness extension of the 
displacement surface on the airfoil surface. We have also carried out solutions 
neglecting the trailing-edge corrections and both the wake curvature and thickness 
terms. This version of the theory, (F) , includes only the displacement effect on the 
airfoil surface and is, therefore, equivalent to the formulation employed in the BGKJ 
(Refs. 12 & 15) method and its derivatives (Refs. 13 Sc 14). For the two 
supercritical cases studied using this latter formulation, we have also obtained 
solutions with the N-C version of our method (labeled (G) in the Table 2B). In 
addition, for the two RAE 2822 cases considered in Tables 2A, B, we have also 
included available results from the RAE non-conservative method (Ref. 17). The RAE 
method accounts for both wake-thickness and wake-curvature but not trailing-edge 
interaction effects. 

The results for the subcritical RAE 2822 case are listed in Table 2A. 
Comparisons of the inviscid (A) and full viscous (B) solutions illustrate the large 
effect of the boundary layer on lift. Even for this subcritical case, the presence 
of the boundary layer causes a nearly 1/3 decrease in lift. The results given in 
lines (C), (D) , and (E) illustrate the effect of the individual terms in the matching 
conditions. Comparisons of the results in lines (A) -(E) indicate that the largest 
effect is caused by the displacement thickness on the airfoil surface with the other 
terms also producing significant effects. The influence of wake thickness on drag 
and lift and of wake curvatures on lift are the most pronounced of these secondary 
influences. The effect of wake thickness on lift was unanticipated since we would 
expect predominately symmetric effects from this source. However, on reflection, 
this behavior is not surprising. The neglect of the wake -thickness term is seen to 
lead to an increase in pressure at the trailing edge which in turn, induces a 
significant increase in boundary- layer thickness near the trailing edge. The 
increase is largest on the upper surface because of the larger initial boundary- layer 
thickness on this surface. It is this differential increase in boundary- layer 
thickness caused by the neglect of the wake-thickness term that leads to the observed 
l^ft reduction. We also call attention to the insensitivity of the momentum drag, 
^^^^ erences i n the theoretical model and to the good agreement between the 
two evaluations of drag in the full viscous solution (B) . The good agreement evident 
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TABLE 2A - RAE 2822 AIRFOIL AT M*, = 0,676 a = 1.06° 
R e = 5.7 x 10 6 X T = 0.11 


Mode 

Theoretical Model 

C D 

C D 

0O 

C L 

CpTE 

® 

Inviscid 

0 

0 

0.571 

0.438 

© 

Full Viscous 

0.0084 

0o0083 

0.431 

0.234 

© 

(b) Less TE 
Corrections 

0.0080 

0.0083 

0,416 

0.234 

® 

(b) Less Wake Curva- 
ture 

0.0081 

0.0083 

'0.454 

0.234 

© 

(b) Less Wake Thick- 
ness 

0.0069 

0.0082 

0o 398 

0.289 


(b) Less TE Correc- 
tions, Wake Curvature 
& Thickness 

0.0064 

0.0082 

0.399 

0,278 

© 

RAE (NC) (Ref 17) 

0.0081 

0.0083 

0,430 

(0»420) 

0,226 


in this case is an encouraging indication of the overall accuracy of the present 
method. The underpredicticn of the integrated drag and the poor agreement with C 

JJ°° 

when the wake- thickness terras are suppressed (E), are indications of the importance 
of these terms. The comparisons in Table 2A indicate that the formulation based on 
airfoil displacement thickness only, underpredicts the drag by 25% and the lift by 7% 
in this case. 

The results of the N-C RAE method of Ref. 17 (C ! ), which is equivalent to the 
formulation employed in case (C) of the present method, are in good agreement. Two 
values of the lift coefficient are given for the RAE method. The higher value 
results from a version of their method that employs considerable numerical smoothing 
of the wake curvature terms near the trailing edge. The smaller value, obtained with 
less smoothing, is in better agreement with the results of case (C) of the present 
method. It was reported in Ref. 17 that convergence difficulties were experienced in 
this latter version and its use was not recommended. This difficulty is likely 
related to the neglect of trailing edge interaction effects in their formulation. 
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A similar series of calculations were carried out for the same airfoil at a 
higher free-stream mach number, resulting in a supercritical flow. The results are 
listed in Table 2B . In this case the integrated drag includes wave drag and so is no 
longer equal to the momentum drag, C . The difference between C _ and C is a 
relatively accurate measure of the wave drag in the full viscous solution (case (B)), 
but not in the other cases because of the inadequacies of the viscous formulation 
clearly evident in the subcritical results (Table 2A) . The neglect of wake and 
trailing edge interaction effects in this supercritical case leads to a 35% 
underprediction of the drag coefficient (case (F)) compared to the full viscous 
solution. We have also repeated case (F) with the N-C formulation (Line (G) in Table 
2B). The N-C scheme leads to a weaker, more forward shock wave but to a higher drag. 
This extra drag arises from a spurious mass generation at the shock wave caused by 
N-C differencing. This spurious drag increase is also evident in the comparison of 
the RAE results (case (C 1 )) with the present, equivalent formulation using F-C 
differencing (case (C)). 

The solutions for the pressure distribution on the airfoil surface are given in 
Fig. 14. Included in the figure are the inviscid solution (A), the full viscous 
solution (B), and the viscous solution less the trailing correction and both wake 
terms (F) . Both F-C and N-C calculations of the latter case are included. These 
results illustrate the very large effect of the boundary layer on the flow field in 
supercritical conditions even at the relatively large Reynolds number of the 
computations. The boundary layer drives the shock wave forward from X = 0.80 to X = 
0.55, significantly reduces the pressure level on the upper surface, and increases it 
on the lower surface of the airfoil. These effects produce a reduction in the lift 
coefficient by nearly a factor of two. The neglect of the wake and trailing-edge 
contributions to the matching conditions is seen to drive the shock wave forward by 
about 5% and to reduce the lift by 10% and drag by 35% below the full viscous 
solution (see Table 2B) . The switch to a N-C scheme is seen to drive the shock wave 
even further forward. 

The details of the pressure distribution near the trailing edge for the full 
viscous solution at the higher Mach number, M^ = 0.725, is plotted on an expanded 
scale in Fig. 15. Included in this figure are the "composite" and "outer" solutions 
computed from the full viscous theory. The composite solution contains contributions 
from the "outer" and "inner" solution as given by Eq. (66) which account for the 
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TABLE 2B - RAE 2822 AIRFOIL AT 0.725, a = 2.3°, 
R e = 6.5 X 10 6 X T = 0.03 


Mode 

Theoretical Model 

C D 

C D 

CO 

C L 

C P 

TE 

® 

Inviscid 

- 

- 

1.300 

0,5265 


Full Viscous 

0.0110 

0.0098 

0.726 

0.222 

© 

(b) Less TE Corrrections 

0.0101 

0.0096 

0.700 

0o 227 

® 

(b) Less Wake Curvature 

0.0112 

0.0100 

0.758 

0.224 

© 

(?) Less Wake Thickness 

0.0080 

0.0093 

0.661 

0.276 

® 

(B) Less TE Corrections, 
Wake Curvature & Thickness 

0.0072 

0.0043 

0.658 

0.271 

© 

(f) With Nonconservative 
Differencing 

0.0081 

0.0092 

0.648 

0.269 

© 

RAE (NC) 5 

0.0119 

- 

0.699 

0.219 


pressure variations across the boundary layer and wake. The outer solution clearly 
exhibits the wake pressure jump imposed as part of the viscous matching conditions. 
The pressure jump is largest at the trailing edge and rapidly decays away from the 
trailing edge both in the wake and on the airfoil surface. The relatively large jump 
in pressure at the trailing edge (AC^ 0.2) is indicative of the importance of 
normal pressure gradients near trailing edges. The crossing of the outer solution on 
the airfoil is a typical property of the viscous wake solution which clearly 
distinguishes it from jet-flap solutions with positive blowing coefficient. 

The RAE 2822 is an airfoil with moderate rear loading. A similar series of 
calculations have also been carried out for a more highly rear- loaded supercritical 
airfoil developed by the NASA-Langley Research Center. The calculations for this 
case were carried out for a free-stream Mach number of M^ = 0.768 and an angle of 
attack of a = -0.151, resulting in a lift coefficient in the full viscous solution of 
C L = 0.852. 

The solutions for this case are listed in Table 3. The airfoil has a nonzero 
trailing-edge thickness which requires an estimate for the base-drag contribution to 
C db - This contribution amounted to 7 counts of drag for all cases considered in 
Table 3 . 
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■ i . ■ A - INVISCID (F-C) 

. B - FULL VISCOUS (F-C) 

— — ^ — — F — FULL VISCOUS (F-C) LESS WAKE AND TE CORRECTIONS 
G - FULL VISCOUS (N-C) LESS WAKE AND TE CORRECTIONS 



Figure 14 RAE 2822 — Effect of the Wake and Trailing Edge Corrections on the Theoretical 
Pressure Distribution at = 0.725, a = 2.3°, R 0 = 6.5 x 10 6 (Xj = 0.03) 


The shock position in the inviscid solution at this incidence and Mach number 
appears to be downstream of the trailing edge and, because of this, the inviscid 
solution could not be computed with the circle plane code used in the present study. 
The viscous solution could be determined without difficulty. The inclusion of the 
boundary layer in the viscous solutions drove the shock wave forward to 73% chord, 
again demonstrating the very large effect of the boundary layer on the flow field at 
transonic speeds. 

The results in Table 3 indicate that the trailing-edge correction and the 
variation of profile drag with the theoretical model are much larger than in previous 
cases. We also note that the neglect of trailing-edge correction and the wake terms 
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Figure 15 RAE 2822 — Details of the Pressure Distribution Near the Trailing Edge at 
= 0 .725, a = 2.3°, R 0 = 6.5 x 10 6 (X y = 0.03) 

(F) leads to a 30% an underprediction of drag and the use of a N-C formulation leads 
to a partially compensating increase of drag of about 10%. 

5.3 THEORETICAL PREDICTION OF PROFILE DRAG AT SUBCRITICAL CONDITIONS 

At low Mach number subcritical conditions, the total drag of an airfoil is equal 
to the profile drag given by the sum of skin friction and pressure drag. As 
discussed previously, the profile drag can be evaluated either by integrating the 
surface pressure and shear-stress distributions or from momentum considerations, by 
C = 20 . The total drag can also be estimated from the Squire-Young approximation 
based on momentum thickness at the trailing edge. The difference between the 
integrated drag, C nT} , and the more accurate momentum drag, C_ , is a useful measure 
of the error of the theory. In Table 4 we compare the theoretical results for the 
integrated and momentum drags for several airfoils at a variety of subcritical flow 
conditions. In this table C^g is the total integrated drag, C^ is the integrated 
drag less the base drag, C^ is momentum drag evaluated from C^^ = 20^ and C^gy 
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total drag evaluated from the Squire-Young formula. The results in the table also 
include experimental data based on wake-rake measurements (Refs. 70-73). 


TABLE 3 ADVANCED SUPERCRITICAL AIRFOIL AT M = 0.768 a = -0.151° 


R = 7.7 X 10 6 X_, = 0.28 

e T 


Mode 

Theoretical Model 

C DB 

C D 

OO 

C L 

CpTE 

© 

Inviscid 

NOT AVAILABLE 

® 

Full viscous 

0.0122 

0.0082 

0.852 

0.137 

© 

( 5 ) Less TE corrections 

0.0097 

0.0078 

0.808 

0.156 

® 

( 5 ) Less wake curvature 

0.0131 

0.0091 

0.871 

0.120 

© 

( 5 ) Less wake thickness 

0.0105 

0.0079 

0.835 

0.155 

© 

® Less TE corrections, 
wake curvature & 
thickness 

0.0088 

0.0079 

0.821 

0.168 

© 

(f) With nonconservative 
differencing 

0.0098 

0.0074 

0.768 

0.155 


The agreement between the integrated and momentum drags for the RAE 2822, GK1 

and the NASA LaRC airfoils is strong evidence that the present method produces 

accurate prediction of the absolute level of drag at subcritical conditions. The 

difference between C n and for the NASA airfoil arises from a thrust on the finite 

base of this airfoil. The agreement between the integrated and momentum drags is 

noticeably poorer for the NACA 0012 airfoil. This airfoil has a significant base 

(Ay/c = 0.0026) and a relatively large trailing-edge angle of 16^. These geometric 

features make the solution for the integrated drag very sensitive to the pressure 

levels near the trailing edge. The under-prediction of the integrated drag in this 

case is very likely caused by neglect of strong-interaction effects associated with 

the trailing-edge angle. The ratios of trailing-edge thickness to boundary- layer 

momentum thickness (Ay/9 ) is 1.2 for the NASA airfoil and 0.4 for the NACA 0012 

lf« 

cases in Table 4, indicating that the approximations used for the base pressure is 
adequate for these cases. The theoretical results for the Squire-Young evaluation of 
the drag is nearly identical to the drag based on far-field momentum thickness in all 
cases considered in the table. The momentum drag, C n , which is thought to be the 

Jjoo 
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more accurate theoretical prediction, is also in very good agreement with the data 
from the RAE, Langely, and Ohio State University wind tunnels. The drag from the NAE 
tunnel seems to be consistently high by about 10 counts compared to the present 
theory. The good agreement observed in other cases (compare, in particular, the Ohio 
State and NAE data for the NACA 0012 airfoil) suggests this is due to a consistent 
error in the NAE measurements . 


TABLE 4 PROFILE DRAG RESULTS 


Case 

Present Theory 

Experiment 

C D 

C DB 

C Dco 

C °SY 

C D Exp Source 

RAE 2822 ^=0.676 

C T =0.576 R =5.7 X 10 6 X =0.11 
ij e I 

— 

0.0087 

0.0085 

0.0085 

0.0085 

RAE 

(Ref. 

70) 

RAE 2822 M^O.676 

C = -0.121 R =5.7 X 10 6 X =0.11 
L e T 

- 

0.0081 

0.0082 

0.0083 

0.0079 

RAE 

(Ref. 

70) 

GKI M oo =0 .511 

C =0.431 R =21.5 X 10 6 X =0.10 
L e T 

— 

0.0063 

0.0064 

0.0064 

0.0078 

NAE 

(Ref. 

73) 

GKI M^O.622 

C T =0.458 R =21.5 X 10 6 X =0.10 
L e T 

— 

0.0065 

0.0066 

0.0066 

0.0070 

NAE 

(Ref. 

73) 

NASA LRC M = 0.78 

C =0.42 R =7.7 X lo 6 X =0.28 
L e T 

0.0079 

0.0070 

0.0070 

0.0069 

0.0072 

NASA 

LRC 

NACA 0012 M =0.575 

OO 

C T =0.006 R =4.7 X 10 6 X =0.10 
L e 1 

0.0079 

0.0072 

0.0082 

0.0082 

0.0081 

Ohio 

State 

Univ. 

(Ref. 

72) 

NACA 0012 M ra =0.693 

C T =0.017 R =22.1 X 10 6 X =0.05 
L e T 

0.0065 

0.0057 

0.0067 

0.0068 

0.0078 

NAE 

(Ref. 

71) 
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In Table 5 the profile drag prediction of the full viscous theory for the NACA 
0012 airfoil tested at Ohio State University are compared with parallel wake 
solutions obtained with the present computer code and with the predictions of the 
BGKJ and Carlson methods presented in Ref* (72)* As discussed above, the momentum 
drag from the full viscous solution is in good agreement with the experimental data* 
The integrated drag shows a significant error for reasons previously mentioned. The 
theoretical pressure distribution is compared with experimental data in Fig. 16. The 
agreement is good, except near the trailing edge, where the theory indicates a more 
positive pressure. This is consistent with the error in integrated drag. 


TABLE 5 NACA 0012 AT ZERO LIFT 

M=0.575 R =4.68 X 10 6 X =0.10 
e 1 


Theory 

C D 

C DB 

Cn 

•Woo 

Cj 

U SY 

Present - Full viscous 

0.0079 

0.0072 

0.0082 

0.0082 

Present - Parallel wake 

0.0063 

0.0054 

0.0081 

0.0081 

BGKJ (Ref. 72) 

0.0047 

0.0036* 

- 

- 

Carlson (Ref. 72) 

0.0051 

0.0042* 

- 

- 

Experiment (Ref. 72) 

0.0081 


(* Obtained by subtracting AC = 0.0009 from C D ) 

B 


Neglect of the wake-thickness term leads to substantial increase in the error in 
integrated drag but does not affect the solution for momentum drag. A significant 
part of the discrepancy between the two evaluations of drag is due to the 
contribution of the base, which in this case amounts to seven counts of thrust. We 
note the very large underprediction of the drag by both the BGKJ and Carlson methods. 
With the base drag included, these methods are seen to underpredict the drag by a 
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factor of two. The base-drag contribution was not included in the original solutions 
presented in Ref. (72). We have corrected these solutions for base drag using the 
base-drag increments determined from our parallel-wake solution in order to provide a 
uniform basis for comparison. 

The results in this section demonstrated that the absolute value of the profile 
drag can be adequately determined from the integrated drag of the full viscous 
solution. This was demonstrated in Table 4 for airfoils with small trailing edge 
angles. The trailing edge angles of the GK1 . NASA-LaRC RAE 2822, and NACA 0012 
airfoil are 0, 3.87, 8.67, 15.96 degrees, respectively. This differs from the 
conclusion reached by Smith and Cebeci (Ref. 74) in a study of boundary- layer methods 
for predicting profile drag. In their study they concluded that profile drag 
predictions based on integrated shear and pressures were inadequate, and they 
recommend exclusive use of the Squire-Young formula. The present study indicated 
that Smith and Cebeci ? s difficulties with integrated drag were due entirely to the 
neglect of the wake- thickness terms in their boundary- layer model. Calculation of 
the Joukowski profile with the present method showed relatively good agreement 
between the integrated and momentum drags. The results are summarized in Table 6. 

We note that the present results for integrated drag obtained with the 
wake-thickness terms neglected agree well with the results of Ref. (74) and that both 
methods underpredict the drag by about a factor of two for the 30% thick airfoil. 

The results in this section demonstrated the importance of the base drag 
contribution for certain airfoils with open trailing edges. The contribution of the 
base drag was an important element in obtaining favorable agreement between theory 
and experiment of the NASA LaRC airfoil. The poorer agreement obtained the the NACA 
0012 is indicative of the importance of strong interaction effects associated with 
the large trailing-edge angle of this airfoil. Further improvements in the 
prediction of profile drag of airfoils with sizeable trailing edge angles can be 
achieved by incorporating an appropriate strong interaction solution into the 
theoretical formulation. 


75 




Figure 16 NACA 0012 - Pressure Distribution at Zero Lift at = 0.575, 
Re = 4.68 x 10 6 (X T =0.10) 


TABLE 6 PROFILE DRAG FOR JOUKOWSKI AIRFOILS AT ZERO INCIDENCE 
= 0.05 R g = 10 X 10 6 X T = 0.10 

Present Theory Theory of Ref. (73) 


t/c 

C f 

C D 

C °sv 

0.10 

0.0058 

0o0064 

(0.0063) 

0.0066 

0.15 

0.0060 

0.0070 

(0.0067) 

0.0074 

0 o 20 

0.0061 

0.0075 

(0.0069) 

0.0081 

0.25 

0.0062 

0.0082 

(0.0069) 

0.0092 

0.30 

0.0063 

0.0091 

(0.0063) 

0.0105 


c f 

C D 


0.0057 

0.0062 

0.0065 

0.0058 

0.0065 

0.0072 

0.0060 

0.0067 

0.0080 

0.0061 

0.0065 

0.0090 

0.0062 

0.0054 

0.0102 


( ) Denotes results obtained with parallel wake formulation 
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5.4 COMPARISONS OF TRANSONIC RESULTS WITH EXPERIMENTAL DATA 

In this section we present comparisons between theoretical computations and 
wind-tunnel data for the GK1 the NASA LaRC, and the RAE 2822 airfoils. The 
comparisons considered in this section include a number of cases of supercritical 
flow with shock waves. The theoretical method described in this study does not 
provide for a proper asymptotic description of the shock-wave boundary- layer 
interaction process. A rational analysis of shock-wave boundary- layer interactions 
in trubulent flow should account for the penetration of the shock wave into the 
boundary layer and normal pressure gradients across the boundary layer. In the 
present study, normal pressure-gradient effects are only accounted for near the 
trailing edge. At shock waves, the present method reduces to a standard, 
free- interaction type analysis employing a Prandtl boundary- layer description of the 
shear layer. This theoretical scheme does not permit shock waves, with the attendant 
discontinuous pressure distribution, to impinge on the boundary layer. A 
discontinuity in pressure or edge velocity in the present method would lead to a 
delta function behavior in the surface source velocity in the viscous matching 
conditions (since v^ % dU^/ds) and to a breakdown of the present method. The results 
of extensive computations suggest this does not occur in the present method. 
Instead, interaction between the inviscid and boundary- layer regions leads to an 
adjustment in the solution that avoids this type of singular breakdown. Within the 
present method, the behavior near shock wave appears to be as sketched in Fig. 17. 



Figure 17 Schematic of a Shock-Wave Boundary Layer Interaction in the Present 
Formulation 
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In this model, the approach of the shock wave to the boundary layer causes a rapid 
buildup of displacement thickness which, in turn, generates a narrow band of 
compression waves. These compressions act to spread the shock wave and prevent a 
discontinuity from impinging on the boundary layer. The interaction spreads the 
pressure jump over a distance, AX~ 6? / y/ M T -1 , where 6 T and M are the displacement 
thickness and local Mach number just upstream of the shock wave. For the cases 
considered in the present study, the interaction length is typically in the order of 
AX~ 0.005. Since the minimum grid spacing on the fine mesh is about twice this 
value, it is clear that the present computations do not resolve the flow near the 
shock wave. Nevertheless, results presented in this section show that the present 
method does lead to reasonable predictions of shock wave/boundary layer interactions. 
It has been known for some time, from experiments on airfoils, that a turbulent 
boundary layer reduces the pressure rise across a shock wave to about half of that 
predicted by the shock conditions for a normal shock wave. This behavior is 
predicted reasonability well in the present method and, in addition, also gives a 
good overall prediction the pressure distribution near the shock wave. 

GK1 Airfoil . - The first group of results are for the GK1 supercritical airfoil, 
the airfoil is 11.5% thick, has a cusped trailing edge, and a moderate degree of aft 
camber. The data are from tests at the NAE transonic wind tunnel at Ottawa (Ref. 73) 
at a Reynolds number of 21.5 million with the tunnel walls set at 20.5% porosity. 
The airfoil was aerodynamical ly smooth and was tested with natural transition. The 
location of the transition point was not determined in the experiments. The 
calculations were carried out with transition fixed at 10% chord which appears to be 
a reasonable estimate, considering the high Reynolds number of the test. The 
theoretical solutions were found to be insensitive to small changes in the 
transition location. Recent studies (Ref. 75) have indicated the presence of 
significant wall interference in the NAE wind tunnel. The studies indicate the need 
for both downwash and blockage corrections to the angle of attack and free-stream 
Mach number with the corrections varying with both free-system Mach number and lift 
coefficient. Unfortunately, corrections were not yet available for the data 
considered in this section. Therefore, in the present comparison, we avoid the need 
for angle of attack correction by comparing theory and experiment at the same value 
of the lift coefficient. Small Mach number shifts are applied to the data to obtain 
agreement in the shock positions. The quantity, a appearing in the Figure is the 

u 

geometric angle of attack of the airfoil in the experiment. 
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The comparison of the surface pressure distribution for two subcritical and two 
supercritical cases are given in Figs. 18 to 21. We note the good agreement of the 
pressure distribution near the trailing edge with experiment and the very good 
predictions of the pressure levels behind the shock wave. The Mach number shift of 
AM^ = 0.008 applied to the data of Fig. 21 is consistent with levels of the blockage 
correction predicted in Ref. 75. The experimental drags are high compared to the 
theoretical value in all cases with the largest discrepancy arising in the = 0.821 
case. The difference in this case is surprisingly high considering the relatively 
close agreement of the pressure distribution. as discussed previously, the 
differences in drag seem to be due to experimental error. Some evidence, for this was 
presented in Ref. 76. 
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Figure 18 GK 1 — Pressure Distribution at = 0.511, C. =0.431, R =21.5x10® 
(X T = 0.10) e 
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Figure 19 GK1 - Pressure Distribution at = 0.622, C, = 0.442, R n = 21.5 (X x = 0.10) 

L. G I 

The solution for the drag polar is compared with experiment and with solutions 
to the Navier-Strokes equations* in Fig. 22. The theoretical solutions were obtained 
with a small Mach number shift of M = -0.005. The present solutions give a 
reasonable prediction of the overall shape of the drag polar but the level of the 
data is about 15% higher than the computation. The poor agreement of the 
Navier-Strokes solutions with both the data and the present computations is most 
likely due to poor spatial resolution since only about 50 points were employed on the 
airfoil in the calculation. 

NASA-La RC Supercritical Airfoil . - This airfoil is a 10% thick, heavily 
rear-loaded, supercritical airfoil with a trailing-edge angle of 3.87° that was 
designed and tested at the NASA LaRC. The airfoil was tested in the 8 foot transonic 
pressure tunnel at a Reynolds number Re = 7.7 X 10^ and a nominal Mach number of M = 

oo 

0.78. Transition was fixed with roughness strip at 28% chord. 

* Unpublished results supplied by Dr. G. S. Deiwert of NASA Ames Research Center. 
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Figure 22 GK1 - Drag Polar at M = 0.752, R e = 21.5 x 10 6 , (X T = 0.10) 


O 

This tunnel is known to be subject to relatively large tunnel-wall interference 

effects, but, unfortunately, wall corrections have not been determined. We applied a 

blockage correction of AM = -0.012 to improve the overall agreement of shock 

00 

position. 

The calculations for the drag polar are compared with wind-tunnel data in Fig. 
23. Theoretical solutions from the full viscous theory are presented along with 
solutions obtained with the wake and trailing-edge correction terms deleted from the 
matching conditions. In the latter case solutions are presented for both F-C and N-C 
difference schemes. The results of the full viscous theory show good agreement with 
experiment over a wide range of lift coefficients. The results show that neglect of 
the wake and trailing-edge terms leads to noticeable underpredictions of drag and 
that N-C differencing leads to a substantial overprediction of drag when shock waves 
are present. 
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Figure 23 NASA LaRC - Drag Polar at M* = 0.768, R 0 = 7.7 x 10 6 , (Xj = 0.28) 


In Fig. 24 we compare the theoretical and experimental pressure distribution for 
the highest lift coefficient of the drag polar. The comparisons are carried out at 
the same lift coefficient (C = 0.94) because of the uncertainties in the effective 

J_i 

angle of attack in the wind tunnel. Note the relatively good agreement in the 
pressure distribution near the shock wave. Although the boundary- layer 
approximations used in the present theory are certainly not valid near shock-wave 
interaction regions, the present method does produce reasonable solutions in these 
regions. The results also show a smooth pressure distribution near the trailing 
edge. The elimination of all numerical smoothings from the present method has 
eliminated the wiggles that appeared in solutions obtained with earlier versions 
(Refs. 5 & 8) of the method. The slightly higher pressures in the theoretical 
solution between the shock wave and the trailing edge are likely due to the neglect 
of strong interaction effects associated with the nonzero, trailing-edge angle of 
this airfoil. 
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In Fig. 25 we present the theoretical solutions for the skin friction and 
displacement thickness on the upper surface of the airfoil. The discontinuities near 
the leading edge are due to transition from laminar to turbulent flow at 28% chord. 
The skin friction shows a nearly discontinuous drop to zero at the shock wave 
indicating that the flow is close to shock-induced separation and probably close to 
stall. We also call attention to the steep and nearly discontinuous jump in 
displacement thickness through the shock wave and the large increase in displacement 
thickness between the shock wave and trailing edge. 

RAE 2822 Airfoil . - The airfoil is 12% thick supercritical section with a 
moderate amount of rear-loading and a trailing edge angle of 8.67°. The airfoil was 
subject to extensive testing (Ref. 70) in the RAE 8x6 foot transonic tunnel in 
which measurements of both surface pressures and boundary layer development were 
carried out. The examples considered in this report were carried out at Reynolds 
numbers from 2.7 x 10^ to 6.5 x 10^ and at two nominal Mach numbers. M = 0.676 and 

oo 

M^ = 0.73, with the lower Mach number resulting in subcritical flow and the higher 
value yielding supercritical flow. The supercritical cases covered a range of 
incidence resulting in a variation of shock strength from moderate to strong-enough 
to cause shock- induced boundary- layer separation. Transition was fixed at 11% chord 
in the lower Mach number (sub-critical) cases and at 3% chord in the higher Mach 
number (supercritical) cases. The momentum and displacement thicknesses, shape 
factor, and skin friction were determined from velocity-profile measurements made at 
a number of locations on the model, and the drag was determined from wake rake 
measurements. All theoretical solutions were carried out at the measured lift 

coefficients. The supercritical cases employed a small blockage correction of AM = 

0.003 which is consistent with recommendations^ of Ref. 18. 

The pressure distribution for the two subcritical cases are compared with 
experimental data in Figs. 26 and 27. These cases were for a Mach number of M^ = 

0.676, Reynolds number Re = 5.7 x 10^, and = 0.566 and = -0.121. The agreement 
with experimental data is generally excellent over the entire airfoil 


# It was actually suggested in Ref. 18 that AM^ = 0.004 be used as a blockage 
correction. However, this reference became available only after the present 
calculations are completed. The slight difference of AM = 0.001 did not warrant a 

oo 

repeat of the computations . 
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Figure 27 RAE 2822 — Pressure Distribution at = 0.676, C. = -0.121, R = 5.7 x 10 


for both cases, with a tendency of the theory to predict a slightly higher pressure on 
the upper surface near the trailing edge. This trend appears in all the RAE 2822 
cases considered, and is likely due to total pressure losses across the shock wave 
that are not accounted for in the potential flow formulation. As previously 
discussed, the prediction of drag for these two subcritical cases was also in very 
good agreement with experimental data (see Table 4). 

The pressure distribution for three supercritical cases tested at a Reynolds 

number of Re = 6.5 x 10^ and a nominal Mach number of M =0.73 are compared with 

00 

experiment in Figs. 28, 29 and 30. Both F-C and N-C solutions are shown in the 
comparisons. In addition, in the last case for C T = 0.803, we have included a 
solution carried out at the experimentally quoted Mach number to illustrate the 
effect of a small Mach number shift (AM^ = 0.003) on the solution. The F-C solution 
clearly shows the best agreement with experiment. In all three cases the agreement 
between theory and experiment is excellent on the lower surface and the shock-wave 
position and strength on the upper surface are well predicted in the F-C solutions. 
This set of results also shows very good agreement with experiment for the pressure 
rise through the shock wave. Both the theoretical and experimental results indicate 
that the pressure rise across the shock wave is only about one-half of that required 
by the normal shock-wave relations. The overall levels of the pressure distribution 
on the upper surface of the airfoil are also reasonably will predicted. The small 
discrepancies between theory and experiment over the forward part of the upper 
surface is probably due to both the roughness strip used to fix transition and 
wall-interference effects. The overprediction of the pressure on the upper surface 
near the trailing edge in Figs. 28 and 29 is likely due to effects associated with 
the nonzero trailing-edge angle not accounted for in our trailing-edge solution. The 
discrepancy is somewhat larger in the higher lift case of Fig. 30, probably because 
of the stronger shock wave and closeness of the boundary layer to separation. The 
poorer agreement may be due to shortcomings in the lag entrainment method at high 
values of the shape factor. 

Since the flow field in these last three cases is supercritical, the total drag 
is given by the sum of wave drag and profile drag. We have previously shown that the 
profile drag of this airfoil can be well predicted by integration of the pressure and 
shear stress. Since the present method seems to give a good prediction of shock 
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Figure 28 RAE 2822 — Pressure Distribution at = 0.728, C 









strength, it can also be expected to give useful predictions of wave drag and total 
drag. This expectation is borne out in the comparison of the drag polar in Fig. 31. 
The theoretical solutions are integrated drags determined from the surface pressure 
and shear stress distributions. We should point out that useful estimates for the 
wave drag can be determined by subtracting the momentum drag (equal to 20^) from the 
total drag computed by the surface integrations. We note that the F-C solution is in 
good agreement with the measured values of drag. The small underprediction of the 
theory evident in the figure amounts to no more than five counts of drag. About half 
of the difference is due to the use of first-order differencing (c^ = 0) in the 
inviscid solution. The remaining discrepancy can be associated with the slight 
overestimate of the pressures on the upper surface near the trailing edge that was 
pointed out earlier. These results clearly show that the N-C scheme leads to an 
inferior prediction of both shock wave location and drag. 



Figure 31 RAE 2822 - Drag Polar at M m a 0.73, R 0 = 6.5 x 10 6 , (X T = 0.03) 
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The theoretical solution (F-C) for the boundary- layer development on the upper 
surface of the airfoil is compared with experimental measurements for the two 
subcritical cases in Figs. 32 and 33 and for two supercritical cases in Figs. 34 and 
35. Included in these figures are the solutions for the displacement and momentum 
thicknesses (6*, 0), shape factor (h) , skin friction (c^), and, for the two 
supercritical cases, the solution for the source velocity, (V^Cs)) on the airfoil 
surface. 

The two subcritical flows in Figures 32 and 33 involve mild pressure gradients 
and boundary layers not close to separation. The discontinuities near the leading 
edge are due to transition from laminar to turbulent flow. Both theory and 
experiment show rapid increases of displacement thickness and shape factor near the 
trailing edge that are typical of airfoil flows. The agreement between theory and 
experiment for all quantities is good for both cases. The good agreement observed 
just downstream of transition indicates that the transition jump conditions employed 
in the theoretical model were adequate. The close agreement between the 
theoretical and experimental values of the momentum thickness at the trailing edge is 
consistent with the good predictions of profile drag in these two cases. 

The boundary- layer development for the two higher- lift, supercritical cases is 
given in Figs. 34 and 35. The overall agreement between theory and experiment is 
relatively good over the entire upper surface aside from a tendency to slightly 
underpredict 6*, 0, and h between the shock wave and the trailing edge. The higher 
lift case considered in Fig. 35 contains data from two different probes. The 
differences between the two data sets provide some measure of the uncertainty in the 
data. The experimental displacement thickness shows a very rapid and nearly 
discontinuous rise through the shock wave, followed by a more gradual but larger 
growth toward the trailing edge. The theory gives a good representation of this 
behavior in both cases. The solution for the skin friction shows the characteristic 
behavior near shock waves; a rapid, nearly discontinuous fall followed by a gradual 
increase resulting in the development of a minimum in skin friction at the shock 
wave. This behavior is reflected in the solution of the shape factor which shows a 
very steep rise followed by a more gradual decrease behind the shock wave. The 
agreement between theory and experiment for c^ and h is quite good except for the 
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points immediately behind the shock wave in the C = 0.803 case in Fig. 35. The 
shock strength in the case is relatively large (The local Mach number upstream of the 
shock wave was M^ = 1.31), and the boundary layer is on the verge of separating at 
the foot of the shock wave. 

The theoretical solutions for the corrected source velocity for the C T = 0.743 
and = 0.803 cases are shown in Figs. 34d and 35d, respectively. These solutions 
show a large and highly localized increase of the source velocity at the shock wave 
followed by a more gradual, but equally significant, growth toward the trailing edge. 
These results clearly show the highly localized nature of the shock wave/boundary 
layer interaction in the present theoretical model The interacting boundary- layer 
description basically truncates the delta function behavior that would otherwise 
arise at the shock impingement point were interaction not allowed in the formulation. 
The resulting source velocity displays a highly peaked but finite distribution which 
increases in amplitude with increasing shock strength (Compare Figs. 34d and 35d.) 

The final case considered is for M^ = 0.753, = 0.743, Re = 6.2 x 10^, and X^, 

“ 0.03. The main difference from the previous cases is a higher Mach number and a 
much lower Reynolds number. These conditions lead to a more rearward and stronger 
shock wave (M^ = 1.35) which, combined with the lower Reynolds number, leads to 
massive separation at the foot of the shock wave. The pressure distribution and 
boundary- layer parameters are compared with experiment in Figs. 36 and 37, 
respectively. The agreement between the theoretical and experimental pressure 
distributions remains excellent on the lower surface but deteriorates noticeably on 
the upper surface. The overall levels of the pressure in front of the shock wave are 
in good agreement on the upper surface except near the leading edge where the 
discrepancy is much larger. However, in the theoretical solution, the shock wave is 
too far aft, its strength is overpredicted, and the pressure levels behind the shock 
wave are overpredicted. The experimental pressures on the upper surface are much 
lower near the trailing edge, showing clear signs of trailing-edge divergence-a 
classical indicator of boundary- layer separation. Although both the theoretical and 
experimental skin friction distributions clearly show separation, pressure divergence 
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Figure 37 RAE 2822 — Boundary Layer Development at M 00 = 0.753, C. = 0.743, R= 6.2 x 10 



at the tailing edge is not apparent in the theoretical results and this must be 
considered a major shortcoming of the therory. The failure of the theory to 
adequately predict the pressure distribution behind the shock wave is likely due to 
an inadequacy in the lag-entrainment method at high values of the shape factor. The 
comparisons of the displacement and momentum thicknesses given in Fig. 35a show that 
the theory gives a reasonably good description of boundary- layer growth in front of 
and through the shock-pressure rise. However, it is seen that the theory seriously 
underpredicts the boundary- layer growth behind the shock wave. It is this 
underprediction of the boundary- layer thickness behind the shock wave that is most 
likely responsible for the observed discrepancy in shock position. The failure of 
the method behind the shock wave is also evident in the comparison of the shape 
factor in Fig. 37b. The comparison of the skin friction is given in Fig. 37c. The 
values of the skin friction in the reversed flow could not be determined in the 
experiment. However, the presence of separation could be detected and such points 
are indicated by the filled symbols in Fig. 37c. As previously noted, the skin 
friction in the theoretical model is set to zero in regions of separation. Note that 
both theory and experiment predict separation behind the shock wave. 

The solution for the surface source velocity in Fig. 37d shows an extremely 
large increase at the shock wave because of separation. In general, this is a rather 
extreme case for the present theory, which was not intended to deal with regions of 
extensive separated flow. Nevertheless, the present theory does seem to provide a 
reasonably good prediction of the overall features of the flow and a surprisingly 
accurate prediction of the drag. The theoretical drag was = 0.0245 compared to 
the experimental value of C = 0.0242. In view of the significant differences 

UD 

between the theoretical and experimental shock strength and location and 
displacement-thickness distributions, this good prediction of the drag should be 
considered fortuitous. 

Collyer and Lock (Ref. 18) recently carried out similar comparisons between the 
RAE 2822 data and results of their method. Their method is similar in many respects 
to the present full viscous -flow theory in that both include the wake-thickness and 
wake -curvature terms in the matching conditions. The Collyer-Lock method differs 
from the present method in three ways: (1) it does not account for the strong 

interaction at the trailing edge; (2) it requires significant numerical smoothings 
particularly near the trailing edge while the present method uses no smoothings; and 
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(3) it employs a special partially-conservative (P-C) scheme while the present 
method uses a standard fully-conservative (F-C) method. The P-C differencing method 
involves an arbitrary parameter, X, that weighs the relative proportion of 
conservative (X = 1) to non-conservative (X = 0) differencing used in the method. 
The position of the shock wave and the value of the drag in the P-C method seemed to 
be sensitive to the value of X. Good agreement with data was demonstrated in Ref. 18 
but for values of X that were varied for the different cases . The present work 
clearly indicated that the best results were obtained with the fully-conservative 
difference scheme in all cases. 

Comparisons of the RAE 2822 data set with Deiwert*s (Ref. 20) Navier-Stokes (NS) 
method were also recently carried out by Swafford in Ref. 77. The Navier-Stokes 
calculations showed very poor agreement with the RAE 2822 data. The poor performance 
of the NS method in this case was very likely due to the coarse grid employed in the 
computations of Ref. 77 (Fourty points were used around the airfoil in the NS 
computations compared to the 160 points used in the present set of calculations.) 
The coarse grid notwithstanding, the NS method required about ten times greater 
computer time than the present interacting boundary- layer method. 
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6. CONCLUDING REMARKS 


In this report we described important extensions to the usual interacting 
boundary- layer theory for computing the viscous effects on airfoils. Our method was 
based on formal asymptotic expansions of the full Reynolds equations of turbulent 
flow for Re-*<». Formal analysis indicated that both wake-induced effects and normal 
pressure-gradient effects near trailing edges contribute to the lowest-order so- 
lution. Consequently, both must be retained along with standard displacement effects 
on the airfoil in order to obtain a consistent formulation. The main contributions 
of the present work were the determination of the local trailing-edge solution and 
the development of systematic procedures for incorporating wake and trail-edge inter- 
action effects into the theoretical description. The inner solution accounted for 
normal pressure gradient effects across the boundary layer in the vicinity of the 
trailing edge and is valid for airfoils with cusped or nearly cusped trailing edges. 

Results presented in this report for a variety of airfoils indicated that the 
new method gave reasonable predictions of the pressure distribution, shock location, 
forces, and moments on airfoils in transonic flows. In particular, the new method 
seems to give very accurate predictions of the absolute levels of the drag. Compara- 
tive studies have indicated that both the wake and trailing-edge corrections terms 
make an important contribution that must be included in the formulation if accurate 
predictions of section characteristics are to be obtained. Neglect of either contri- 
bution leads to noticeable errors in shock position, pressure distribution, lift, and 
drag. In addition, the present study clearly indicated the superiority of conserva- 
tive finite-difference formulations. The use of non-conservative differencing leads 
to substantial over-predictions of the drag and to shock positions that are too far 
forward. The present method is more complete and is clearly an improvements over 
previous boundary- layer type methods. It agrees well with experiment, requires lit- 
tle computing time, and is, therefore, well suited to aerodynamic design. The method 
requires about 10 minutes on an IBM 370-168 computer to obtain converged solutions 
for unseparated cases. This is about three times that required for a corresponding 
calculation of an inviscid flow. These computer times are for unsmoothed versions of 
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the theory. Faster run times can be achieved with numerical smoothings but only at 
the expense of noticeably degraded solutions. 

The present method does not properly treat the details of shock wave/boundary 
layer interactions. These interactions produce significant normal pressure gradients 
in the boundary layer that are completely ignored in the simple interacting boundary- 
layer equations used in our approach. Nevertheless, the new method seems to give 
adequate predictions of the pressure distributions in the shock-wave interaction 
zones. It is a well known experimental fact that the pressure rise through a shock 
wave is about half that required by the normal shock-wave relations. This behavior 
is well represented by the present method which also seems to give a consistently 
good prediction of the pressure level behind the shock wave. The consequences of 
using a crude description of the shock wave/boundary layer interaction process do not 
appear to be great. The overall increase in displacement thickness through the shock 
wave is reasonably well predicted by the present method, and this seems to be the 
most important ingredient in obtaining good overall predictions of the section char- 
acteristics of airfoils in supercritical flow. 

The present method has two other more important deficiencies that should be men- 
tioned. Namely, it does not properly treat the flow near the trailing edge of an 
airfoil with non-zero trailing edge angles, and it does not adequately deal with 
boundary layer separation. Fortunately, the present method seems to yield reasonably 
good predictions when applied to more general airfoils with non-zero trailing edge 
angles provided the trailing edge angles are small. However, the solution for the 
NACA 0012 airfoil (0,^ = 16°) showed pressures that were somewhat too high near the 
trailing edge and drags that were substantially lower than experimental values. 
These discrepancies were undoubtedly caused by pressure variations across the bounda- 
ry layer in the vicinity of the trailing edge that are not accounted for in our local 
inner solution. These extra pressure variations across the boundary layer are gen- 
erated by the curved streamlines associated with the local wedge flow near the trail- 
ing edge. These effects can be incorporated into the method through an appropriate 
generalization of our inner trailing-edge solution. This improvement of the theory 
would very likely result in better predictions of drag for this type of airfoil and 
should therefore be pursued. Although the present method functions for separated 
flow has produced reasonably good agreement with experiment in some cases, it is 
clearly not satisfactory in this regard. The lag-entrainment method was developed 
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from a data base for attached boundary layers and was not intended to deal with sep- 
arated flows or with boundary layers with large shape factors. Treatment of 
separation in the present method can be improved through improvements in the 
lag-entrainment method (Some effort in this direction has been reported (Ref. 78).) 
or through the use of other integral methods (Ref. 79) specifically designed for sep- 
arated turbulent boundary- layers . In addition, the iteration scheme used to solve 
the coupled inviscid and boundary layer equations loses some effectiveness when sepa- 
ration is present. In these situations, solutions can be obtained only through the 
use of drastically reduced relaxation factors, resulting in very long computer times. 
Improved convergence can probably be achieved through the use of an iteration scheme 
developed by Carter (Ref. 80) and LeBalleur (Refs. 81 & 82) for separated flows. 

The most important result of the present work is the demonstration that a simple 
interacting boundary- layer approach can be effective tool for predicting viscous ef- 
fects on airfoils. Unfortunately, because of uncertainties in the experimental data, 
due primarily to wind-tunnel wall interference, this must be regarded as a tentative 
conclusion. A more clear-cut validation of the new theory would require improved 
wind-tunnel tests with either small wall effects or accurate and well-documented cor- 
rections . 
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7. EPILOGUE 


The present report was completed in early 1980 and the computations reported 
herein were carried out with a preliminary version of the GRUMFOIL computer code des- 
ignated MCMJ-4. Since that time, the GRUMFOIL code has undergone extensive modifica- 
tion to incorporate a variety of improvements made to the theoretical formulation. 
First, the hybrid, SLOR-direct solver, used to solve the full potential equation, was 
replaced with a multi-grid method (Ref. 3) which accelerated convergence by a factor 
of five in computing time. Second, an improved, second order artificial viscosity 
and far field treatment was incorporated into the finite difference scheme used for 
the inviscid flow. These significantly reduced the spatial truncation error, im- 
proved the reliability of the second-order versions of the code, and reduced by 
one-half the number of mesh points required to achieve a given accuracy. Finally, 
the turbulence closure relations employed in the original Green s lag entrainment 
method modified to improve the treatment of separated flow, and the direct method 
employed for the global, viscid-inviscid iteration was replaced with Carter s 
semi-inverse method to permit solution of the resulting turbulent boundary layer 
equations when separation is present. Together, these changes led to a more accu- 
rate, robust, and faster code which is capable of treating flows with extensive 
regions of separated flow and which runs about ten times faster than the original 
MCMJ-4 version. Further details of the new code, designated MCMJ-9N, are given in 
Ref. 83. The user manual in the companion volume, Ref. 69, is based on the new 
MCMJ-9N version of the code. 
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